firemodels · drjfloyd · Dec 20, 2024 · Dec 20, 2024 · Dec 20, 2024
diff --git 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
@@ -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
@@ -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<CIRCLE_BBOX_AREA) THEN ! more efficient to descretize the rectangle
-   ! make area elements nominally square
-   IF ((X2-X1)<(Y2-Y1)) THEN
-      NX = 50
-      BBDX = (X2-X1)/REAL(NX,EB)
-      NY = NINT((Y2-Y1)/BBDX)
-      BBDY = (Y2-Y1)/REAL(NY,EB)
-   ELSE
-      NY = 50
-      BBDY = (Y2-Y1)/REAL(NY,EB)
-      NX = NINT((X2-X1)/BBDY)
-      BBDX = (X2-X1)/REAL(NX,EB)
-   ENDIF
-   DELTA_AREA = BBDX*BBDY
-   DO JJ=1,NY
-      DO II=1,NX
-         XX = X1 + BBDX*(II-0.5_EB)
-         YY = Y1 + BBDY*(JJ-0.5_EB)
-         IF ( ((XX-X0)**2+(YY-Y0)**2)<RAD**2 ) CIRCLE_CELL_INTERSECTION_AREA = CIRCLE_CELL_INTERSECTION_AREA + DELTA_AREA
-      ENDDO
-   ENDDO
-ELSE BBOX_IF ! more efficient to descretize the circle
-   NX = 50
-   NY = 50
-   CX1 = X0-RAD
-   CY1 = Y0-RAD
-   BBDX = 2._EB*RAD/REAL(NX,EB)
-   BBDY = BBDX
-   DELTA_AREA = BBDX*BBDY
-   DO JJ=1,NY
-      DO II=1,NX
-         XX = CX1 + BBDX*(II-0.5_EB)
-         YY = CY1 + BBDY*(JJ-0.5_EB)
-         IF ( ((XX-X0)**2+(YY-Y0)**2)<RAD**2 .AND. &
-                XX>=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
@@ -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
@@ -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