diff --git a/Source/dump.f90 b/Source/dump.f90 index 7488a531a91..d88f39f4077 100644 --- a/Source/dump.f90 +++ b/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). diff --git a/Source/geom.f90 b/Source/geom.f90 index e3a066a4029..22e6cea471c 100644 --- a/Source/geom.f90 +++ b/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,9 +8877,6 @@ 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) @@ -8845,6 +8884,8 @@ SUBROUTINE INIT_CFACE_CELL(NM,ICF,IFACE,CFACE_INDEX,SURF_INDEX,STAGE_FLG,IS_INB, 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' diff --git a/Source/hvac.f90 b/Source/hvac.f90 index b77475d2fed..10949e00e72 100644 --- a/Source/hvac.f90 +++ b/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 diff --git a/Source/init.f90 b/Source/init.f90 index 40997c7313b..6f7a9e58447 100644 --- a/Source/init.f90 +++ b/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 diff --git a/Source/read.f90 b/Source/read.f90 index 9a5d009fe9f..e8e71110377 100644 --- a/Source/read.f90 +++ b/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 diff --git a/Source/type.f90 b/Source/type.f90 index 73a5ea1c4e9..3d4619e0e1a 100644 --- a/Source/type.f90 +++ b/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