FireX: Merge with firemodels master

Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex

@@ -1245,7 +1245,7 @@ \section{2-D Vortex Simulation (\texorpdfstring{\ct{vort2d}}{vort2d})}
 \textbf{Max Gould, NIST SURF student}\\
 \textbf{Ragini Acharya, United Technologies Research Center}\\
 
-\noindent In this section we present another case that demonstrates the second-order accuracy of the FDS transport algorithm. We consider the analytically stable flow field consisting of a single vortex advected by a uniform flow, a test case developed by CERFACS\footnote{Centre Europ\'een de Recherche et de Formation Avanc\'ee en Calcul Scientifique} \cite{cerfacs_test}. Maintaining the geometry of the vortex over time provides a good measure of the order of accuracy of the transport scheme.
+\noindent In this section another case is presented case that demonstrates the second-order accuracy of the FDS transport algorithm. Consider the analytically stable flow field consisting of a single vortex advected by a uniform flow, a test case developed by CERFACS\footnote{Centre Europ\'een de Recherche et de Formation Avanc\'ee en Calcul Scientifique} \cite{cerfacs_test}. Maintaining the geometry of the vortex over time provides a good measure of the order of accuracy of the transport scheme.
 
 \begin{figure}[h!]
    \centering
@@ -1287,33 +1287,32 @@ \section{2-D Vortex Simulation (\texorpdfstring{\ct{vort2d}}{vort2d})}
 u (x,z) &\equiv& U_{0} + \frac{\partial}{\partial z} \Psi_{0} = U_{0} - \frac{\Gamma \ z}{R_{c}^{2}} \exp \left[ - \frac{x^{2} + z^{2}}{2 \ R_{c}^{2}} \right], \\
 w (x,z) &\equiv& - \frac{\partial}{\partial x} \Psi_{0} = \frac{\Gamma \ x}{R_{c}^{2}} \exp \left[ - \frac{x^{2} + z^{2}}{2 \ R_{c}^{2}} \right],
 \end{eqnarray}
-where $u$ and $w$ refer to velocity in the $x$ and $z$-directions, respectively. For our purposes we need only analyze one component of the velocity field. We will focus our attention on the $u$-component of velocity.
+where $u$ and $w$ refer to velocity in the $x$ and $z$-directions, respectively. It is sufficient to analyze the $u$-component of velocity.
 
-We define the computational domain as a two-dimensional square region, $L=0.3112$~m on a side, with periodic boundary conditions. The domain is discretized for a range of square, two-dimensional meshes of $40^{2}$, $80^{2}$, $160^{2}$, and $320^{2}$ grid cells. For the purposes of this test, we set the flow parameters as
+The computational domain is a two-dimensional square region, $L=0.3112$~m on a side, with periodic boundary conditions. The domain is discretized for a range of square, two-dimensional meshes of $40^{2}$, $80^{2}$, $160^{2}$, and $320^{2}$ grid cells. The flow parameters are
 \[
 \begin{array}{lll}
 U_{0} &&= 35 \ \mathrm{m/s} \\
 R_{c} &= L / 20 &= 0.01556 \ \mathrm{m} \\
 \Gamma &= 0.04 \ U_{0} \ R_{c} \ \sqrt{e} &= 0.0359157
 \end{array}
 \]
-The constant flow field and periodic boundary conditions cause the vortex to repeatedly pass through the computational domain. The ``pass-through'' time, $t_{f}$, is defined as the time period required for the stable vortex to return to its original position,
+The constant flow field and periodic boundary conditions cause the vortex to repeatedly pass through the computational domain. The ``pass-through'' time, $t_{\rm f}$, is defined as the time period required for the stable vortex to return to its original position,
 \begin{equation*}
-t_{f} = L / U_{0} \simeq 8.8914 \times 10^{-3} \ \mathrm{s}.
+t_{\rm f} = L / U_{0} \simeq 8.8914 \times 10^{-3} \ \mathrm{s}
 \end{equation*}
-To ensure that the numerical solution converges to the analytical solution, we set the time step, $\dt$, so that the Courant-Friedrichs-Lewy (CFL) number is 0.5.
-
+To ensure that the numerical solution converges to the analytical solution, the time step, $\dt$, is chosen so that the Courant-Friedrichs-Lewy (CFL) number is 0.5.
 
 \begin{figure}[h]
    \begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r}
       \includegraphics[height=2.2in]{SCRIPT_FIGURES/vort2d_80_uzgraph} &
       \includegraphics[height=2.2in]{SCRIPT_FIGURES/vort2d_160_uzgraph}
    \end{tabular*}
-   \caption[Velocity in the \ct{vort2d} test case]{$u$-velocity along the $z$-axis at $x=0$ plotted for each of the vortex's first four loops through the computational domain. (Left) $80^{2}$ grid cell model. (Right) $160^{2}$ grid cell model.}
+   \caption[Velocity in the \ct{vort2d} test case]{$u$-velocity along the $z$-axis at $x=0$ plotted for each of the vortex's first three passes through the computational domain. (Left) $80^{2}$ grid cell model. (Right) $160^{2}$ grid cell model.}
    \label{fig_vort2d_axisvelocity}
 \end{figure}
 
-A plot of $u$-velocity values just along the $z$-axis provides a simple characterization of the vortex geometry. The extent to which this geometry changes over time provides a qualitative measure of the accuracy of the transport algorithm. Figure \ref{fig_vort2d_axisvelocity} displays these plots for two different grid resolutions. Each line represents a plot taken for a different number of flow-through times such that the red lines represent the vortex after it has undergone the most passes through the computational domain while the green lines represent the vortex in the initial phase. The broken black line represents the analytical solution. As the vortex undergoes more passes through the computational domain, its velocity profile diverges further and further from the analytical profile. While divergence still occurs on the finer mesh, the extent to which it diverges after the same number of flow-through times is significantly smaller.
+A plot of $u$-velocity along the $z$-axis provides a simple characterization of the vortex geometry. The extent to which this geometry changes over time provides a qualitative measure of the accuracy of the transport algorithm. Figure~\ref{fig_vort2d_axisvelocity} displays these plots for two different grid resolutions. Each curve represents the state of the vortex at each flow-through time. The black curve represents the analytical solution for which the numerical solution has been initialized. As the vortex undergoes more passes through the computational domain, its velocity profile diverges from the analytical profile. While divergence still occurs on the finer mesh, the extent to which it diverges after the same number of flow-through times is significantly smaller.
 
 \begin{figure}[h!]
    \centering
@@ -1322,7 +1321,7 @@ \section{2-D Vortex Simulation (\texorpdfstring{\ct{vort2d}}{vort2d})}
     \label{fig_vort2d_error}
 \end{figure}
 
-The rate with which the simulated profiles converge to the analytical one defines the order of accuracy of the numerical scheme. In Fig.~\ref{fig_vort2d_error} we plot the rms error of the numerical solution as a function of the grid resolution. The three colored curves represent the rms error at three different pass-through times. The broken and solid black lines represent the plot gradient corresponding to first and second order error respectively. While the error increases with each flow-through time, the gradients of the lines are roughly parallel to the solid black line, indicating second-order accuracy of the numerical scheme.
+The rate with which the simulated profiles converge to the analytical one defines the order of accuracy of the numerical scheme. In Fig.~\ref{fig_vort2d_error} the error of the numerical solution is plotted as a function of the grid resolution. The three colored curves represent the error at three different pass-through times. The broken and solid black lines represent the plot gradient corresponding to first and second order error respectively. While the error increases with each flow-through time, the gradients of the lines are roughly parallel to the solid black line, indicating second-order accuracy of the numerical scheme.
 
 To analyze the stability of the vortex at times other than discrete multiples of the pass-through time, consider the time-dependent potential field and its corresponding $u$-velocity component:
 \begin{equation}
@@ -1333,7 +1332,7 @@ \section{2-D Vortex Simulation (\texorpdfstring{\ct{vort2d}}{vort2d})}
 \label{eqn_uvel_timedep}
 u (x,z,t) = U_{0} - \Psi_{0} \ \frac{z}{R_{c}^{2}} \exp \left[\frac{2 \ U_{0} \ x \ t - U_{0}^{2} \ t^{2}}{2 \ R_{c}^{2}} \right].
 \end{equation}
-In Fig.~\ref{fig_vort2d_pointvelocity}, we show the $u$-velocity at a single point, on the lower left fringe of the vortex, for two different mesh resolutions, $80^{2}$ and $160^{2}$. A mesh resolution of $320^{2}$ grid cells produces a plot (not shown) that is, to the eye, a perfect match out to four pass-through times.
+In Fig.~\ref{fig_vort2d_pointvelocity}, the $u$-velocity at a single point on the lower left fringe of the vortex is plotted for two different mesh resolutions, $80^{2}$ and $160^{2}$. A mesh resolution of $320^{2}$ grid cells produces a plot (not shown) that is, to the eye, a perfect match out to four pass-through times.
 
 \begin{figure}[h!]
    \begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r}
@@ -1491,7 +1490,7 @@ \subsection{Fire Plume using Constant Specific Heat Ratio (\texorpdfstring{\ct{f
 \subsection{Evaporation with Constant Specific Heat Ratio (\texorpdfstring{\ct{water_evap_1_const_gamma}}{water\_evap\_1\_const\_gamma})}
 \label{water_evap_1_const_gamma}
 
-This test case is a replica of \ct{water_evaporation_1} in Sec.~\ref{water_evaporation_1} using constant specific heat ratio.  A 1 \si{m^3} box is initially filled with dry air at 200 \si{\degreeCelsius} and mono-disperse water droplets totaling 0.01 kg in mass initially at 20 \si{\degreeCelsius}.  In this case, all the water evaporates.  The final gas temperature may be computed from energy conservation.  Pressure change may be computed from the ideal gas law.  Results are shown in Fig.~\ref{water_evap_1_const_gamma_plots}.  Details of the expected results may be found in \ct{water_evap_1_const_gamma.m} in the \ct{Utilities/Matlab/scripts/} directory in the FDS repository \cite{FDS-SMV_repository}.
+This test case is a replica of \ct{water_evaporation_1} in Sec.~\ref{water_evaporation_1} using constant specific heat ratio.  A 1 \si{m^3} box is initially filled with dry air at 200 \si{\degreeCelsius} and mono-disperse water droplets totaling 0.01 kg in mass initially at 20 \si{\degreeCelsius}.  In this case, all the water evaporates.  The final gas temperature may be computed from energy conservation.  Pressure change may be computed from the ideal gas law.  Results are shown in Fig.~\ref{water_evap_1_const_gamma_plots}.  Details of the expected results may be found in \ct{water_evap_1_const_gamma.py} in the \ct{Utilities/Python/scripts/} directory in the FDS repository \cite{FDS-SMV_repository}.
 
 \begin{figure}[ht!]
 \noindent
@@ -7828,6 +7827,7 @@ \chapter{Outputs}
 \section{Statistical Quantities}
 
 \subsection{RMS, Co-Variance, and Cross-Correlation (\texorpdfstring{\ct{rms\_cov\_corr}}{rms\_cov\_corr})}
+\label{rms_cov_corr}
 
 FDS can output the root mean square (RMS), co-variance, and cross-correlation for both point and line \ct{DEVC} outputs.  To test these outputs a 1 m$^3$ box with open sides and a 10 cm grid size is defined with two inlet vents centered on adjacent faces.  This results in two orthogonal flow streams that collide at the center of the box and exit diagonally.  Within the diagonal portion of the flow, are placed point measurements for the FDS outputs of the u-velocity RMS, the u-velocity/w-velocity co-variance, and the u-velocity/w-velocity cross-correlation.
 

Source/ccib.f90

@@ -107,7 +107,7 @@ MODULE CC_SCALARS
           CC_CUTCELL_VELOCITY,CC_CUTFACE_VELOCITY,CC_RESTORE_UVW_UNLINKED,&
           CHECK_CFLVN_LINKED_CELLS,ADD_Q_DOT_CUTCELLS,CFACE_THERMAL_GASVARS,&
           CFACE_PREDICT_NORMAL_VELOCITY,COMPUTE_LINKED_CUTFACE_BAROCLINIC,&
-          COPY_CC_UNKH_TO_HS, COPY_CC_HS_TO_UNKH, COPY_UNST_DM_TO_CART, CUTFACE_VELOCITIES, &
+          COPY_CC_UNKH_TO_HS, COPY_CC_MUNKH_TO_UNKH, COPY_UNST_DM_TO_CART, CUTFACE_VELOCITIES, &
           GET_CFACE_OPEN_BC_COEF,GET_FN_DIVERGENCE_CUTCELL,GET_OPENBC_TANGENTIAL_CUTFACE_VEL,&
           GET_CUTCELL_DDDT,GET_H_CUTFACES,GET_H_MATRIX_CC,GET_H_GUARD_CUTCELL,GET_CRTCFCC_INT_STENCILS,GET_RCFACES_H, &
           GET_CC_MATRIXGRAPH_H,GET_CC_IROW,GET_CC_UNKH,GET_CUTCELL_HP, GET_LINKED_FV, GET_PRES_CFACE_BCS, &
@@ -21750,44 +21750,57 @@ SUBROUTINE NUMBER_UNKH_CUTCELLS(FLAG12,NM,IPZ,NUNKH_LC)
 RETURN
 END SUBROUTINE NUMBER_UNKH_CUTCELLS
 
-! ------------------- COPY_CC_HS_TO_UNKH ----------------------------
+! ------------------- COPY_CC_MUNKH_TO_UNKH ----------------------------
 
-SUBROUTINE COPY_CC_HS_TO_UNKH(NM)
-
-INTEGER, INTENT(IN) :: NM
+SUBROUTINE COPY_CC_MUNKH_TO_UNKH
 
 ! Local Variables:
-INTEGER :: NOM,ICC,II,JJ,KK,IOR,IW,IIO,JJO,KKO
+INTEGER :: NOM,ICC,IW,IIO,JJO,KKO,II,JJ,KK
 TYPE (OMESH_TYPE), POINTER :: OM
-
+TYPE (BOUNDARY_COORD_TYPE), POINTER :: BC
 ! Loop over external wall cells:
 EXTERNAL_WALL_LOOP: DO IW=1,N_EXTERNAL_WALL_CELLS
 
    WC=>WALL(IW)
    EWC=>EXTERNAL_WALL(IW)
    IF (WC%BOUNDARY_TYPE/=INTERPOLATED_BOUNDARY) CYCLE EXTERNAL_WALL_LOOP
 
-   BC => BOUNDARY_COORD(WC%BC_INDEX)
-   II = BC%II
-   JJ = BC%JJ
-   KK = BC%KK
-   IOR = BC%IOR
    NOM = EWC%NOM
    OM => OMESH(NOM)
 
    ! This assumes all meshes at the same level of refinement:
    KKO=EWC%KKO_MIN
    JJO=EWC%JJO_MIN
    IIO=EWC%IIO_MIN
-
+   BC => BOUNDARY_COORD(WC%BC_INDEX)
+   II = BC%II; JJ = BC%JJ; KK = BC%KK;
    ICC=CCVAR(II,JJ,KK,CC_IDCC)
-
    IF (ICC > 0) THEN ! Cut-cells on this guard-cell Cartesian cell.
-      MESHES(NM)%CUT_CELL(ICC)%UNKH(1) = INT(OM%HS(IIO,JJO,KKO))
+      CUT_CELL(ICC)%UNKH(1) = OM%MUNKH(IIO,JJO,KKO)
    ELSE
-      MESHES(NM)%CCVAR(II,JJ,KK,CC_UNKH) = INT(OM%HS(IIO,JJO,KKO))
+      CCVAR(II,JJ,KK,CC_UNKH) = OM%MUNKH(IIO,JJO,KKO)
    ENDIF
 
+   ! ! Loop over all cells in mesh NOM that correspond to this boundary face (supports grid refinement):
+   ! DO KKO = EWC%KKO_MIN, EWC%KKO_MAX
+   !    DO JJO = EWC%JJO_MIN, EWC%JJO_MAX
+   !       DO IIO = EWC%IIO_MIN, EWC%IIO_MAX
+
+   !          ! Check if this cell in mesh NOM is a cut-cell:
+   !          ICC = MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_IDCC)
+
+   !          IF (ICC > 0) THEN
+   !             ! Copy to cut-cell storage in mesh NOM:
+   !             MESHES(NOM)%CUT_CELL(ICC)%UNKH(1) = OM%MUNKH(IIO,JJO,KKO)
+   !          ELSE
+   !             ! Copy to regular cell storage in mesh NOM:
+   !             MESHES(NOM)%CCVAR(IIO,JJO,KKO,CC_UNKH) = OM%MUNKH(IIO,JJO,KKO)
+   !          ENDIF
+
+   !       ENDDO
+   !    ENDDO
+   ! ENDDO
+
 ENDDO EXTERNAL_WALL_LOOP
 
 ! Loop over external wall cells:
@@ -21801,7 +21814,7 @@ SUBROUTINE COPY_CC_HS_TO_UNKH(NM)
 ENDDO EXTERNAL_WALL_LOOP2
 
 RETURN
-END SUBROUTINE COPY_CC_HS_TO_UNKH
+END SUBROUTINE COPY_CC_MUNKH_TO_UNKH
 
 ! ------------------- COPY_CC_UNKH_TO_HS ----------------------------
 

Source/main.f90

@@ -1743,15 +1743,19 @@ SUBROUTINE GLOBAL_MATRIX_REASSIGN(FORCE_REASSIGN)
             CALL ULMAT_SOLVER_SETUP(NM)
          ENDDO
          CALL STOP_CHECK(1)
-      CASE (UGLMAT_FLAG,GLMAT_FLAG)
-         IF(ALLOCATED(ZONE_SOLVE)) CALL FINISH_GLMAT_SOLVER
-         CALL GLMAT_SOLVER_SETUP(1)
-         CALL STOP_CHECK(1)
-         CALL MESH_EXCHANGE(3) ! Exchange guard cell info for CCVAR(I,J,K,CGSC) -> HS.
-         CALL GLMAT_SOLVER_SETUP(2)
-         CALL MESH_EXCHANGE(3) ! Exchange guard cell info for CCVAR(I,J,K,UNKH) -> HS.
-         CALL GLMAT_SOLVER_SETUP(3)
-         CALL STOP_CHECK(1)
+   CASE (UGLMAT_FLAG,GLMAT_FLAG)
+      IF(ALLOCATED(ZONE_SOLVE)) CALL FINISH_GLMAT_SOLVER
+      CALL GLMAT_SOLVER_SETUP(-1) ! Initialize EWC_TYPE, copy wall types to HS
+      CALL STOP_CHECK(1)
+      CALL MESH_EXCHANGE(3)       ! Exchange guard cell info for IS_WALLT -> HS
+      CALL GLMAT_SOLVER_SETUP(0)  ! Process coarse faces, copy updated wall types to HS
+      CALL MESH_EXCHANGE(3)       ! Re-exchange guard cell info for IS_WALLT -> HS
+      CALL GLMAT_SOLVER_SETUP(1)
+      CALL MESH_EXCHANGE(3)       ! Exchange guard cell info for CCVAR(I,J,K,CGSC) -> HS
+      CALL GLMAT_SOLVER_SETUP(2)
+      CALL MESH_EXCHANGE(3)       ! Exchange guard cell info for CCVAR(I,J,K,UNKH) -> HS
+      CALL GLMAT_SOLVER_SETUP(3)
+      CALL STOP_CHECK(1)
    END SELECT
    OBST_CREATED_OR_REMOVED = .FALSE.
 ENDIF