diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index 33cb624fe2f..5e626ea8f4c 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -1739,7 +1739,7 @@ \subsubsection{Circular Vents} \begin{lstlisting} &SURF ID='BURNER', MASS_FLUX(1)=0.02, SPEC_ID(1)='PROPANE', TAU_MF(1)=0.01 / &VENT XB=-0.6,0.6,-0.6,0.6,0,0, XYZ=0,0,0, RADIUS=0.5, SURF_ID='BURNER', - SPREAD_RATE=0.05 / + SPREAD_RATE=0.05/ \end{lstlisting} The \ct{XB} coordinates designate the orientation of the vent. In this case, the extent of the area specified by \ct{XB} is large enough to contain the entire circle. Note also in this example that the parameter \ct{SPREAD_RATE} causes the fire to spread outward at a rate of 0.05~m/s. The mass flux of propane through the vent is plotted in Fig.~\ref{circ_burn}. Notice that the mass flux increases following a ``t-squared'' profile. This is what is expected of a fire which spreads radially at a linear rate. In this case, the fire reaches the \ct{RADIUS} of the circle in 10~s, as expected. Note also that the parameter \ct{TAU_MF} indicates that the fuel should ramp up quickly once the flame front reaches a given grid cell. In other words, \ct{TAU_MF} controls the local ramp-up of fuel; the \ct{SPREAD_RATE} controls the global ramp-up. Following the ramp-up, the fuel flows at a rate equal to the area of the circle times the mass flux of fuel per unit area. Even if the circle is crudely resolved on a coarse grid, the fuel flow rate will be adjusted to produce the desired value governed by the circular vent. \begin{figure}[h] @@ -8700,9 +8700,9 @@ \subsection{Freezing the Output Value, Example Case: \ct{hrr_freeze}} \subsection{Example Case: Heat Release Rate of a Spreading Fire} \label{spreading_fire} -In this example, a fire spreads radially from a single point as directed by the parameters \ct{SPREAD_RATE} and \ct{XYZ} on a \ct{VENT} line. Usually, the user specifies the heat release rate per unit area (\ct{HRRPUA}) for each burning surface cell on the corresponding \ct{SURF} line, but in this case, a specific time \ct{RAMP} for the {\em total} heat release rate is specified. The following input lines show how the user-specified \ct{RAMP} called \ct{'HRR'} controls the total HRR of the growing fire. The key point is that the user-specified {\em total} HRR is divided by the area of burning surface, and this heat release rate per unit area is imposed on all burning cells. Regardless of the fact that the spreading fire reaches a barrier and is stopped from spreading radially, the user-specified \ct{RAMP} controls the HRR, as shown in Fig.~\ref{spreading_fire_hrr}. +In this example, a fire spreads radially from a single point as directed by the parameters \ct{SPREAD_RATE} and \ct{XYZ} on a \ct{VENT} line. Usually, the user specifies the heat release rate per unit area (\ct{HRRPUA}) for each burning surface cell on the corresponding \ct{SURF} line, but in this case, a specific time \ct{RAMP} for the {\em total} heat release rate is specified. The following input lines show how the user-specified \ct{RAMP} called \ct{'HRR'} controls the total HRR of the growing fire. The key point is that the user-specified {\em total} HRR is divided by the area of burning surface, and this heat release rate per unit area is imposed on all burning cells. Normally FDS will adjust a mass flux input (\ct{MASS_FLUX}, \ct{HRRPUA} ,etc.) input to account for any differences in the area of the \ct{VENT} as specified with \ct{XB} and the area is it is actually resolved on the grid. In this case we are using control functions to determine the heat release rate. In this case the control logic is directly computing the required flux based on the area as resolved so no additional correction is needed. When false, the \ct{AREA_ADJUST} flag prevents any additional adjustment. Regardless of the fact that the spreading fire reaches a barrier and is stopped from spreading radially, the user-specified \ct{RAMP} controls the HRR, as shown in Fig.~\ref{spreading_fire_hrr}. \begin{lstlisting} -&VENT SURF_ID='FIRE', XB=1.0,7.0,1.0,7.0,0.0,0.0, XYZ=3.0,3.0,0.0, SPREAD_RATE=0.1 / +&VENT SURF_ID='FIRE', XB=1.0,7.0,1.0,7.0,0.0,0.0, XYZ=3.0,3.0,0.0, SPREAD_RATE=0.1, AREA_ADJUST=F / &SURF ID='FIRE', HRRPUA=1., RAMP_Q='HRRPUA RAMP' / &RAMP ID='HRRPUA RAMP', T=0, F=0, CTRL_ID_DEP='HRRPUA CTRL' / @@ -13621,6 +13621,7 @@ \section{\texorpdfstring{{\tt VENT}}{VENT} (Vent Parameters)} \hline Keyword & Type & Description & Units & Default \\ \hline \hline \endhead +\ct{AREA_ADJUST} & Logical & Section~\ref{spreading_fire} & & \ct{F} \\ \hline \ct{COLOR } & Character & Section~\ref{info:colors} & & \\ \hline \ct{CTRL_ID } & Character & Section~\ref{info:activate_deactivate} & & \\ \hline \ct{DB } & Character & Section~\ref{info:VENT_Basics} & & \\ \hline diff --git a/Source/cons.f90 b/Source/cons.f90 index 53b0a48dd25..62254ceceec 100644 --- a/Source/cons.f90 +++ b/Source/cons.f90 @@ -378,7 +378,7 @@ MODULE GLOBAL_CONSTANTS TYPE (MPI_COMM), ALLOCATABLE, DIMENSION(:) :: MPI_COMM_CLOSE_NEIGHBORS !< MPI communicator for the a given mesh and its neighbors INTEGER, ALLOCATABLE, DIMENSION(:) :: MPI_COMM_NEIGHBORS_ROOT !< The rank of the given mesh within the MPI communicator INTEGER, ALLOCATABLE, DIMENSION(:) :: MPI_COMM_CLOSE_NEIGHBORS_ROOT !< The rank of the given mesh within the MPI communicator -INTEGER, ALLOCATABLE, DIMENSION(:) :: COUNTS,DISPLS,COUNTS_10,DISPLS_10,COUNTS_20,DISPLS_20 +INTEGER, ALLOCATABLE, DIMENSION(:) :: COUNTS,DISPLS,COUNTS_10,DISPLS_10,COUNTS_20,DISPLS_20,COUNTS_VENT,DISPLS_VENT ! Time parameters @@ -753,6 +753,9 @@ MODULE GLOBAL_CONSTANTS REAL(EB), ALLOCATABLE, DIMENSION(:) :: EXTERNAL_RAMP LOGICAL, ALLOCATABLE, DIMENSION(:) :: EXTERNAL_CTRL +! VENT array +REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: VENT_TOTAL_AREA + END MODULE GLOBAL_CONSTANTS diff --git a/Source/func.f90 b/Source/func.f90 index 31d6331e61d..bfe47c4e2e9 100644 --- a/Source/func.f90 +++ b/Source/func.f90 @@ -3248,65 +3248,122 @@ SUBROUTINE BLOCK_MESH_INTERSECTION_VOLUME(NM,X1,X2,Y1,Y2,Z1,Z2,VOLUME_ADJUST) END SUBROUTINE BLOCK_MESH_INTERSECTION_VOLUME -!> \brief Estimate area of intersection of circle and rectangle +!> \brief Estimate area of intersection of circle and grid cell !> \param X0 x-coordinate of the center of the circle (m) !> \param Y0 y-coordinate of the center of the circle (m) !> \param RAD Radius of the circle (m) -!> \param X1 Lower x-coordinate of the rectangle (m) -!> \param X2 Upper x-coordinate of the rectangle (m) -!> \param Y1 Lower y-coordinate of the rectangle (m) -!> \param Y2 Upper y-coordinate of the rectangle (m) +!> \param X1 Lower x-coordinate of the grid cell (m) +!> \param X2 Upper x-coordinate of the grid cell (m) +!> \param Y1 Lower y-coordinate of the grid cell (m) +!> \param Y2 Upper y-coordinate of the grid cell (m) REAL(EB) FUNCTION CIRCLE_CELL_INTERSECTION_AREA(X0,Y0,RAD,X1,X2,Y1,Y2) REAL(EB), INTENT(IN) :: X0,Y0,RAD,X1,X2,Y1,Y2 -INTEGER :: NX,NY,II,JJ -REAL(EB) :: DELTA_AREA,XX,YY,CIRCLE_BBOX_AREA,RECT_BBOX_AREA,CX1,CY1,BBDX,BBDY +INTEGER :: FC +REAL(EB) :: XC,YC,XP1,XP2,YP1,YP2 +REAL(EB) :: RC,R2,THETA -CIRCLE_BBOX_AREA = 4._EB*RAD*RAD -RECT_BBOX_AREA = (X2-X1)*(Y2-Y1) +R2 = RAD**2 +FC = 0 -CIRCLE_CELL_INTERSECTION_AREA = 0._EB +! No overlap -BBOX_IF: IF (RECT_BBOX_AREA=X1 .AND. XX<=X2 .AND. & - YY>=Y1 .AND. YY<=Y2 ) CIRCLE_CELL_INTERSECTION_AREA = CIRCLE_CELL_INTERSECTION_AREA + DELTA_AREA - ENDDO - ENDDO -ENDIF BBOX_IF +IF ((X2 < X0-RAD) .OR. (X1 > X0+RAD) .OR. (Y2 < Y0-RAD) .OR. (Y1 > Y0+RAD)) THEN + CIRCLE_CELL_INTERSECTION_AREA = 0._EB + RETURN +ENDIF + +! Count corners inside circle +RC = (X1-X0)**2+(Y1-Y0)**2 +IF (RC <= R2) FC=IBSET(FC,0) +RC = (X1-X0)**2+(Y2-Y0)**2 +IF (RC <= R2) FC=IBSET(FC,1) +RC = (X2-X0)**2+(Y2-Y0)**2 +IF (RC <= R2) FC=IBSET(FC,2) +RC = (X2-X0)**2+(Y1-Y0)**2 +IF (RC <= R2) FC=IBSET(FC,3) + +SELECT CASE(FC) + ! Grid cell surrounds the circle: area of circle + CASE (0) + CIRCLE_CELL_INTERSECTION_AREA = PI*R2 + ! One corner in the circle: chord + area of triangle + CASE(1,2,4,8) + SELECT CASE(FC) + CASE(1) + XC=X1 + YC=Y1 + YP1=SQRT(R2-(X1-X0)**2)+Y0 + XP1=SQRT(R2-(Y1-Y0)**2)+X0 + CASE(2) + XC=X1 + YC=Y2 + YP1=Y0-SQRT(R2-(X1-X0)**2) + XP1=SQRT(R2-(Y2-Y0)**2)+X0 + CASE(4) + XC=X2 + YC=Y2 + YP1=Y0-SQRT(R2-(X2-X0)**2) + XP1=X0-SQRT(R2-(Y2-Y0)**2) + CASE(8) + XC=X2 + YC=Y1 + YP1=SQRT(R2-(X2-X0)**2)+Y0 + XP1=X0-SQRT(R2-(Y1-Y0)**2) + END SELECT + THETA = 2._EB*ASIN(0.5_EB*SQRT((XC-XP1)**2+(YC-YP1)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + 0.5_EB*ABS(XC-XP1)*ABS(YC-YP1) + ! Two corners in the circle: chord + trapezoid + CASE(3) + XP1=SQRT(R2-(Y1-Y0)**2)+X0 + XP2=SQRT(R2-(Y2-Y0)**2)+X0 + THETA = 2._EB*ASIN(0.5_EB*SQRT((XP1-XP2)**2+(Y2-Y1)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + 0.5_EB*(Y2-Y1)*(XP1+XP2-2._EB*X1) + CASE(6) + YP1=Y0-SQRT(R2-(X1-X0)**2) + YP2=Y0-SQRT(R2-(X2-X0)**2) + THETA = 2._EB*ASIN(0.5_EB*SQRT((X2-X1)**2+(YP2-YP1)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + 0.5_EB*(X2-X1)*(2._EB*Y2-YP1-YP2) + CASE(9) + YP1=SQRT(R2-(X1-X0)**2)+Y0 + YP2=SQRT(R2-(X2-X0)**2)+Y0 + THETA = 2._EB*ASIN(0.5_EB*SQRT((X2-X1)**2+(YP2-YP1)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + 0.5_EB*(X2-X1)*(YP1+YP2-2._EB*Y1) + CASE(12) + XP1=X0-SQRT(R2-(Y1-Y0)**2) + XP2=X0-SQRT(R2-(Y2-Y0)**2) + THETA = 2._EB*ASIN(0.5_EB*SQRT((XP1-XP2)**2+(Y2-Y1)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + 0.5_EB*(Y2-Y1)*(2._EB*X2-XP1-XP2) + ! Three corners in the circle: chord + irregular pentagon (two rectangles and a triangle) + CASE(7) + YP1=Y0-SQRT(R2-(X2-X0)**2) + XP1=SQRT(R2-(Y1-Y0)**2)+X0 + THETA = 2._EB*ASIN(0.5_EB*SQRT((XP1-X2)**2+(YP1-Y1)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + (Y2-Y1)*(XP1-X1) + (X2-XP1)*(Y2-YP1) + & + 0.5_EB*(X2-XP1)*(YP1-Y1) + CASE(11) + YP1=SQRT(R2-(X2-X0)**2)+Y0 + XP1=SQRT(R2-(Y2-Y0)**2)+X0 + THETA = 2._EB*ASIN(0.5_EB*SQRT((XP1-X2)**2+(YP1-Y2)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + (Y2-Y1)*(XP1-X1) + (X2-XP1)*(YP1-Y1) + & + 0.5_EB*(X2-XP1)*(Y2-YP1) + CASE(13) + YP1=SQRT(R2-(X1-X0)**2)+Y0 + XP1=X0-SQRT(R2-(Y2-Y0)**2) + THETA = 2._EB*ASIN(0.5_EB*SQRT((XP1-X1)**2+(YP1-Y2)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + (Y2-Y1)*(X2-XP1) + (XP1-X1)*(YP1-Y1) + & + 0.5_EB*(XP1-X1)*(Y2-YP1) + CASE(14) + YP1=Y0-SQRT(R2-(X1-X0)**2) + XP1=X0-SQRT(R2-(Y1-Y0)**2) + THETA = 2._EB*ASIN(0.5_EB*SQRT((XP1-X1)**2+(YP1-Y1)**2)/RAD) + CIRCLE_CELL_INTERSECTION_AREA = 0.5_EB*R2*(THETA-SIN(THETA)) + (Y2-Y1)*(X2-XP1) + (XP1-X1)*(Y2-YP1) + & + 0.5_EB*(XP1-X1)*(YP1-Y1) + ! Entire grid cell is in circle: area of grid cell + CASE(15) + CIRCLE_CELL_INTERSECTION_AREA = (X2-X1)*(Y2-Y1) +END SELECT END FUNCTION CIRCLE_CELL_INTERSECTION_AREA diff --git a/Source/init.f90 b/Source/init.f90 index bc26bd0ebd6..657f1a2668d 100644 --- a/Source/init.f90 +++ b/Source/init.f90 @@ -3331,7 +3331,15 @@ SUBROUTINE INIT_WALL_CELL(NM,I,J,K,OBST_INDEX,IW,IOR,SURF_INDEX,IERR,TT) VT => M%VENTS(WC%VENT_INDEX) - B1%AREA_ADJUST = SF%AREA_MULTIPLIER * VT%INPUT_AREA/VT%FDS_AREA + IF (VT%AREA_ADJUST) THEN + IF (VT%RADIUS > 0._EB) THEN + B1%AREA_ADJUST = SF%AREA_MULTIPLIER * VT%UNDIVIDED_INPUT_AREA/VT%TOTAL_FDS_AREA + ELSE + B1%AREA_ADJUST = SF%AREA_MULTIPLIER * VT%INPUT_AREA/VT%FDS_AREA + ENDIF + ELSE + B1%AREA_ADJUST = SF%AREA_MULTIPLIER + ENDIF IF (B1%AREA_ADJUST<=TWO_EPSILON_EB) B1%AREA_ADJUST = 1._EB IF (VT%CTRL_INDEX > 0) THEN diff --git a/Source/main.f90 b/Source/main.f90 index 191f1814fb6..44459ac1161 100644 --- a/Source/main.f90 +++ b/Source/main.f90 @@ -72,7 +72,7 @@ PROGRAM FDS CHARACTER(MPI_MAX_PROCESSOR_NAME) :: PNAME LOGICAL, ALLOCATABLE, DIMENSION(:) :: LOGICAL_BUFFER_EXTERNAL REAL(EB), ALLOCATABLE, DIMENSION(:) :: REAL_BUFFER_DUCT,REAL_BUFFER_EXTERNAL -REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: REAL_BUFFER_10,REAL_BUFFER_20 +REAL(EB), ALLOCATABLE, DIMENSION(:,:) :: REAL_BUFFER_10,REAL_BUFFER_20,REAL_BUFFER_VENT ! Initialize OpenMP @@ -149,6 +149,10 @@ PROGRAM FDS CALL MPI_INITIALIZATION_CHORES(2) +! Exchange VENT data + +IF (N_VENT_TOTAL > 0) CALL EXCHANGE_VENT_AREA + ! Initialize global parameters CALL INITIALIZE_GLOBAL_VARIABLES @@ -1053,6 +1057,7 @@ SUBROUTINE MPI_INITIALIZATION_CHORES(TASK_NUMBER) ALLOCATE(REAL_BUFFER_DUCT((2+N_TRACKED_SPECIES)*N_DUCTNODES+N_DUCTS)) ALLOCATE(REAL_BUFFER_10(10,NMESHES)) ALLOCATE(REAL_BUFFER_20(20,NMESHES)) + ALLOCATE(REAL_BUFFER_VENT(N_VENT_TOTAL,NMESHES)) IF (READ_EXTERNAL) THEN ALLOCATE(REAL_BUFFER_EXTERNAL(N_RAMP)) @@ -1063,11 +1068,13 @@ SUBROUTINE MPI_INITIALIZATION_CHORES(TASK_NUMBER) ALLOCATE(COUNTS_10(0:N_MPI_PROCESSES-1)) ALLOCATE(COUNTS_20(0:N_MPI_PROCESSES-1)) ALLOCATE(COUNTS_TP(0:N_MPI_PROCESSES-1)) + ALLOCATE(COUNTS_VENT(0:N_MPI_PROCESSES-1)) ALLOCATE(DISPLS(0:N_MPI_PROCESSES-1)) ALLOCATE(DISPLS_10(0:N_MPI_PROCESSES-1)) ALLOCATE(DISPLS_20(0:N_MPI_PROCESSES-1)) ALLOCATE(DISPLS_TP(0:N_MPI_PROCESSES-1)) + ALLOCATE(DISPLS_VENT(0:N_MPI_PROCESSES-1)) ALLOCATE(I_OFFSET(NMESHES)) COUNTS = 0 @@ -1092,6 +1099,8 @@ SUBROUTINE MPI_INITIALIZATION_CHORES(TASK_NUMBER) DISPLS(N) = COUNTS(N-1) + DISPLS(N-1) DISPLS_TP(N) = COUNTS_TP(N-1) + DISPLS_TP(N-1) ENDDO + COUNTS_VENT = COUNTS*N_VENT_TOTAL + DISPLS_VENT = DISPLS*N_VENT_TOTAL COUNTS_10 = COUNTS*10 DISPLS_10 = DISPLS*10 COUNTS_20 = COUNTS*20 @@ -4260,6 +4269,28 @@ SUBROUTINE CHECK_FREEZE_VELOCITY_STATUS END SUBROUTINE CHECK_FREEZE_VELOCITY_STATUS +!> \brief Exchange FDS VENT areas so each process knows the total FDS area + +SUBROUTINE EXCHANGE_VENT_AREA +INTEGER:: NM,NV + +IF (N_MPI_PROCESSES>1) THEN + REAL_BUFFER_VENT = VENT_TOTAL_AREA + CALL MPI_ALLGATHERV(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,REAL_BUFFER_VENT,& + COUNTS_VENT,DISPLS_VENT,MPI_DOUBLE_PRECISION,MPI_COMM_WORLD,IERR) + VENT_TOTAL_AREA = REAL_BUFFER_VENT +ENDIF + +DO NM=1,NMESHES + IF (PROCESS(NM)/=MY_RANK) CYCLE + DO NV=1,MESHES(NM)%N_VENT + MESHES(NM)%VENTS(NV)%TOTAL_FDS_AREA=SUM(VENT_TOTAL_AREA(MESHES(NM)%VENTS(NV)%TOTAL_INDEX,:)) + ENDDO +ENDDO + +END SUBROUTINE EXCHANGE_VENT_AREA + + !> \brief Gather revision dates, etc, from subroutines !> \param REVISION String containing the revision number !> \param REVISION_DATE String containing the date of the last code revision diff --git a/Source/read.f90 b/Source/read.f90 index 729972deb4b..486c1d3374c 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -11699,21 +11699,21 @@ SUBROUTINE READ_VENT USE MATH_FUNCTIONS, ONLY: GET_RAMP_INDEX USE MISC_FUNCTIONS, ONLY: PROCESS_MESH_NEIGHBORHOOD -INTEGER :: N,N_TOTAL,N_EXPLICIT,NM,NNN,IOR,I1,I2,J1,J2,K1,K2,RGB(3),N_EDDY,II,JJ,KK,OBST_INDEX,N_IMPLICIT_VENTS,I_MODE,& +INTEGER :: N,N_TOTAL,NM,NNN,IOR,I1,I2,J1,J2,K1,K2,RGB(3),N_EDDY,II,JJ,KK,OBST_INDEX,N_EXPLICIT,N_IMPLICIT_VENTS,I_MODE,& N_ORIGINAL_VENTS,IC0,IC1 REAL(EB) :: SPREAD_RATE,TRANSPARENCY,XYZ(3),TMP_EXTERIOR,DYNAMIC_PRESSURE,XB_USER(6),XB_MESH(6), & REYNOLDS_STRESS(3,3),L_EDDY,VEL,VEL_RMS,L_EDDY_IJ(3,3),UVW(3),RADIUS 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,GEOM,SOLID_FOUND +LOGICAL :: REJECT_VENT,OUTLINE,GEOM,SOLID_FOUND,AREA_ADJUST 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/ COLOR,CTRL_ID,DB,DEVC_ID,DYNAMIC_PRESSURE,FYI,GEOM,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 @@ -11738,6 +11738,8 @@ SUBROUTINE READ_VENT IF (I_MODE==2) THEN ALLOCATE(MESHES(NM)%VENTS(N_VENT),STAT=IZERO) ; CALL ChkMemErr('READ','VENTS',IZERO) ; VENTS=>MESHES(NM)%VENTS + IF (.NOT. ALLOCATED(VENT_TOTAL_AREA)) ALLOCATE(VENT_TOTAL_AREA(N_VENT_TOTAL,NMESHES),STAT=IZERO) ! Only once per MPI process + CALL ChkMemErr('READ','VENT_TOTAL_AREA',IZERO); VENT_TOTAL_AREA = 0._EB IF (MY_RANK==0 .AND. .NOT.ALLOCATED(ORIGINAL_VENTS)) ALLOCATE(ORIGINAL_VENTS(N_ORIGINAL_VENTS)) ENDIF @@ -11968,6 +11970,8 @@ SUBROUTINE READ_VENT VT=>VENTS(N_VENT) + VT%AREA_ADJUST = AREA_ADJUST + VT%TOTAL_INDEX = N_VENT_TOTAL ! Set basic VENT coordinates VT%I1 = I1 @@ -11992,7 +11996,6 @@ SUBROUTINE READ_VENT VT%Z2_ORIG = XB_USER(6) ! Vent area - IF (ABS(XB_USER(1)-XB_USER(2))<=SPACING(XB_USER(2))) & VT%UNDIVIDED_INPUT_AREA = (XB_USER(4)-XB_USER(3))*(XB_USER(6)-XB_USER(5)) IF (ABS(XB_USER(3)-XB_USER(4))<=SPACING(XB_USER(4))) & @@ -12000,7 +12003,15 @@ SUBROUTINE READ_VENT IF (ABS(XB_USER(5)-XB_USER(6))<=SPACING(XB_USER(6))) & VT%UNDIVIDED_INPUT_AREA = (XB_USER(2)-XB_USER(1))*(XB_USER(4)-XB_USER(3)) - IF (RADIUS>0._EB) VT%UNDIVIDED_INPUT_AREA = PI*RADIUS**2 + ! For circular vent account for VENT XB may only partially overlap circle + IF (RADIUS>0._EB) THEN + IF (ABS(XB_USER(1)-XB_USER(2))<=SPACING(XB_USER(2))) VT%UNDIVIDED_INPUT_AREA = & + CIRCLE_CELL_INTERSECTION_AREA(XYZ(2),XYZ(3),RADIUS,XB_USER(3),XB_USER(4),XB_USER(5),XB_USER(6)) + IF (ABS(XB_USER(3)-XB_USER(4))<=SPACING(XB_USER(4))) VT%UNDIVIDED_INPUT_AREA = & + CIRCLE_CELL_INTERSECTION_AREA(XYZ(1),XYZ(3),RADIUS,XB_USER(1),XB_USER(2),XB_USER(5),XB_USER(6)) + IF (ABS(XB_USER(5)-XB_USER(6))<=SPACING(XB_USER(6))) VT%UNDIVIDED_INPUT_AREA = & + CIRCLE_CELL_INTERSECTION_AREA(XYZ(1),XYZ(2),RADIUS,XB_USER(1),XB_USER(2),XB_USER(3),XB_USER(4)) + ENDIF IF (ABS(VT%X2-VT%X1)<=SPACING(VT%X2) ) VT%INPUT_AREA = (VT%Y2-VT%Y1)*(VT%Z2-VT%Z1) IF (ABS(VT%Y2-VT%Y1)<=SPACING(VT%Y2) ) VT%INPUT_AREA = (VT%X2-VT%X1)*(VT%Z2-VT%Z1) @@ -12318,8 +12329,6 @@ SUBROUTINE READ_VENT ENDIF VT%FDS_AREA = 0._EB - IF (VT%RADIUS>0._EB) VT%INPUT_AREA = 0._EB - I1 = VT%I1 I2 = VT%I2 J1 = VT%J1 @@ -12329,37 +12338,34 @@ SUBROUTINE READ_VENT SELECT CASE(ABS(VT%IOR)) CASE(1) + IF ( VT%RADIUS>0._EB) & + VT%INPUT_AREA = CIRCLE_CELL_INTERSECTION_AREA(VT%Y0,VT%Z0,VT%RADIUS,VT%Y1,VT%Y2,VT%Z1,VT%Z2) DO K=K1+1,K2 DO J=J1+1,J2 - IF ( VT%RADIUS>0._EB) THEN - VT%INPUT_AREA = VT%INPUT_AREA + CIRCLE_CELL_INTERSECTION_AREA(VT%Y0,VT%Z0,VT%RADIUS,Y(J-1),Y(J),Z(K-1),Z(K)) - IF (((YC(J)-VT%Y0)**2+(ZC(K)-VT%Z0)**2)>VT%RADIUS**2) CYCLE - ENDIF + IF (VT%RADIUS>0._EB .AND. ((YC(J)-VT%Y0)**2+(ZC(K)-VT%Z0)**2)>VT%RADIUS**2) CYCLE VT%FDS_AREA = VT%FDS_AREA + DY(J)*DZ(K) ENDDO ENDDO CASE(2) + IF ( VT%RADIUS>0._EB) & + VT%INPUT_AREA = CIRCLE_CELL_INTERSECTION_AREA(VT%X0,VT%Z0,VT%RADIUS,VT%X1,VT%X2,VT%Z1,VT%Z2) DO K=K1+1,K2 DO I=I1+1,I2 - IF ( VT%RADIUS>0._EB) THEN - VT%INPUT_AREA = VT%INPUT_AREA + CIRCLE_CELL_INTERSECTION_AREA(VT%X0,VT%Z0,VT%RADIUS,X(I-1),X(I),Z(K-1),Z(K)) - IF (((XC(I)-VT%X0)**2+(ZC(K)-VT%Z0)**2)>VT%RADIUS**2) CYCLE - ENDIF + IF ( VT%RADIUS>0._EB .AND. ((XC(I)-VT%X0)**2+(ZC(K)-VT%Z0)**2)>VT%RADIUS**2) CYCLE VT%FDS_AREA = VT%FDS_AREA + DX(I)*DZ(K) ENDDO ENDDO CASE(3) + IF ( VT%RADIUS>0._EB) & + VT%INPUT_AREA = CIRCLE_CELL_INTERSECTION_AREA(VT%X0,VT%Y0,VT%RADIUS,VT%X1,VT%X2,VT%Y1,VT%Y2) DO J=J1+1,J2 DO I=I1+1,I2 - IF ( VT%RADIUS>0._EB) THEN - VT%INPUT_AREA = VT%INPUT_AREA + CIRCLE_CELL_INTERSECTION_AREA(VT%X0,VT%Y0,VT%RADIUS,X(I-1),X(I),Y(J-1),Y(J)) - IF (((XC(I)-VT%X0)**2+(YC(J)-VT%Y0)**2)>VT%RADIUS**2) CYCLE - ENDIF + IF ( VT%RADIUS>0._EB .AND. ((XC(I)-VT%X0)**2+(YC(J)-VT%Y0)**2)>VT%RADIUS**2) CYCLE VT%FDS_AREA = VT%FDS_AREA + DX(I)*DY(J) ENDDO ENDDO END SELECT - + VENT_TOTAL_AREA(VT%TOTAL_INDEX,NM) = VT%FDS_AREA ENDDO VENT_LOOP_3 ENDDO MESH_LOOP_2 @@ -12369,6 +12375,7 @@ SUBROUTINE READ_VENT SUBROUTINE SET_VENT_DEFAULTS +AREA_ADJUST = .TRUE. COLOR = 'null' CTRL_ID = 'null' DB = 'null' diff --git a/Source/type.f90 b/Source/type.f90 index a637973273a..0f30e0d1852 100644 --- a/Source/type.f90 +++ b/Source/type.f90 @@ -1480,15 +1480,15 @@ MODULE TYPES TYPE VENTS_TYPE INTEGER :: I1=-1,I2=-1,J1=-1,J2=-1,K1=-1,K2=-1,BOUNDARY_TYPE=0,IOR=0,SURF_INDEX=0,DEVC_INDEX=-1,CTRL_INDEX=-1, & COLOR_INDICATOR=99,TYPE_INDICATOR=0,ORDINAL=0,PRESSURE_RAMP_INDEX=0,TMP_EXTERIOR_RAMP_INDEX=0,NODE_INDEX=-1, & - OBST_INDEX=0 + OBST_INDEX=0,TOTAL_INDEX=-1 INTEGER, DIMENSION(3) :: RGB=-1 REAL(EB) :: TRANSPARENCY = 1._EB REAL(EB), DIMENSION(3) :: TEXTURE=0._EB - REAL(EB) :: X1=0._EB,X2=0._EB,Y1=0._EB,Y2=0._EB,Z1=0._EB,Z2=0._EB,FDS_AREA=0._EB, & + REAL(EB) :: X1=0._EB,X2=0._EB,Y1=0._EB,Y2=0._EB,Z1=0._EB,Z2=0._EB,FDS_AREA=0._EB,TOTAL_FDS_AREA=0._EB,& 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.,GEOM=.FALSE. + LOGICAL :: ACTIVATED=.TRUE.,GEOM=.FALSE.,AREA_ADJUST=.TRUE. CHARACTER(LABEL_LENGTH) :: DEVC_ID='null',CTRL_ID='null',ID='null' ! turbulent inflow (experimental) INTEGER :: N_EDDY=0 diff --git a/Verification/Controls/spreading_fire.fds b/Verification/Controls/spreading_fire.fds index b7b90bf21c8..5f2820a2433 100644 --- a/Verification/Controls/spreading_fire.fds +++ b/Verification/Controls/spreading_fire.fds @@ -14,7 +14,7 @@ &VENT MB='YMIN', SURF_ID='OPEN' / &VENT MB='YMAX', SURF_ID='OPEN' / &VENT MB='ZMAX', SURF_ID='OPEN' / -&VENT SURF_ID='FIRE', XB=1.0,7.0,1.0,7.0,0.0,0.0, XYZ=3.0,3.0,0.0, SPREAD_RATE=0.1 / +&VENT SURF_ID='FIRE', XB=1.0,7.0,1.0,7.0,0.0,0.0, XYZ=3.0,3.0,0.0, SPREAD_RATE=0.1,AREA_ADJUST=F / &SURF ID='FIRE', HRRPUA=1., RAMP_Q='HRRPUA RAMP' / &RAMP ID='HRRPUA RAMP', T=0, F=0, CTRL_ID_DEP='HRRPUA CTRL' /