FireX: Merge with firemodels/master

.gitignore

@@ -45,6 +45,7 @@ Manuals/*/*Guide.pdf
 *.sf
 *.bf
 *.s3d
+*.s3dd
 *.q
 *.sz
 *.smv

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -991,12 +991,11 @@ \subsection{Basics}
 \subsection{Two-Dimensional and Axially-Symmetric Calculations}
 \label{info:2D}
 
-The governing equations solved in FDS are written in terms of a three dimensional Cartesian coordinate system. However, a two dimensional Cartesian or two dimensional cylindrical (axially-symmetric) calculation can be performed by setting the \ct{J} in the \ct{IJK} triplet to 1 on the \ct{MESH} line. For axial symmetry, add \ct{CYLINDRICAL=T} to the \ct{MESH} line, and the coordinate $x$ is then interpreted as the radial coordinate $r$. If more than one mesh is used, all the meshes must be specified as 2-D or \ct{CYLINDRICAL}---you cannot mix 2-D, 3-D and cylindrical geometries. No boundary conditions should be set at the planes $y=\hbox{\tt YMIN=XB(3)}$ or $y=\hbox{\tt YMAX=XB(4)}$, nor at $r=\hbox{\tt XMIN=XB(1)}$ in an axially-symmetric calculation if $r=\hbox{XB(1)=0}$ (Note that \ct{XB(1)} does not have to be 0). For better visualizations, the difference between \ct{XB(4)} and \ct{XB(3)} should be small so that the Smokeview rendering appears to be in 2-D. An example of an axially-symmetric helium plume is given in Sec.~\ref{baroclinic_torque}.
-
-When processing results for a \ct{CYLINDRICAL} simulation, note that integrated output quantities with the \ct{SPATIAL_STATISTIC} attribute apply only to the specified 2-D or cylindrical coordinates. Thus, the cylindrical coordinates define a cylindrical sector, like a slice of cake, even though Smokeview will not render it this way. The fully integrated quantity can be calculated by multiplying the reported value by $2 \pi \, \delta\theta$, where $\delta\theta$ is the difference between \ct{YMAX} and \ct{YMIN} in radians. The values chosen for \ct{YMAX} and \ct{YMIN} do not matter as long as the rendering in Smokeview is to your liking.
+The governing equations solved in FDS are written in terms of a three-dimensional Cartesian coordinate system. However, a two-dimensional Cartesian or two-dimensional cylindrical (axially-symmetric) calculation can be performed by setting the \ct{J} in the \ct{IJK} triplet to 1 on the \ct{MESH} line. For axial symmetry, add \ct{CYLINDRICAL=T} to the \ct{MESH} line, and the coordinate $x$ is then interpreted as the radial coordinate $r$. If more than one mesh is used, all the meshes must be specified as 2-D or \ct{CYLINDRICAL}---you cannot mix 2-D, 3-D and cylindrical geometries. No boundary conditions should be set at the planes $y=\hbox{\tt YMIN=XB(3)}$ or $y=\hbox{\tt YMAX=XB(4)}$, nor at $r=\hbox{\tt XMIN=XB(1)}$ in an axially-symmetric calculation if $r=\hbox{XB(1)=0}$ (Note that \ct{XB(1)} does not have to be 0). For better visualizations, the difference between \ct{XB(4)} and \ct{XB(3)} should be small so that the Smokeview rendering appears to be in 2-D. An example of an axially-symmetric helium plume is given in Sec.~\ref{baroclinic_torque}.
 
 When performing solid phase heat transfer while using a 2-D \ct{CYLINDRICAL} coordinate system, you must designate \ct{GEOMETRY='CYLINDRICAL'} on a surface (\ct{SURF} line) that is facing radially outward (positive $r$ direction) or \ct{GEOMETRY='INNER CYLINDRICAL'} on a surface that is facing radially inward (negative $r$ direction). In the latter instance, you must also specify the \ct{INNER_RADIUS} (m) of the cylinder. For the outer cylindrical boundary, specify an \ct{INNER_RADIUS} if appropriate. Its default value is 0~m. Because your inward and outward facing boundaries might occur at various radii, you must create separate \ct{SURF} lines for each with the appropriate values of \ct{GEOMETRY} and \ct{INNER_RADIUS}. For an obstruction (\ct{OBST}), use \ct{SURF_ID6} to assign individual \ct{SURF ID}s to each of the six faces. Because this is a 2-D simulation, the third and fourth entries representing the \ct{SURF ID}s in the $y$ or angular direction can just be designated \ct{'INERT'}.
 
+When processing results for a \ct{CYLINDRICAL} simulation, note that integrated output quantities with the \ct{SPATIAL_STATISTIC} attribute refer to the volume or surface area of the entire cylinder, not just the wedge. Smokeview renders the wedge as a 2-D slice. The values chosen for \ct{YMAX} and \ct{YMIN} do not matter as long as the rendering in Smokeview is to your liking. For a 2-D, non-cylindrical geometry, spatially integrated quanties shall be output in units of the quantity per unit meter.
 
 \subsection{Multiple Meshes}
 \label{info:multimesh}
@@ -9140,7 +9139,7 @@ \section{Flux Limiters}
 Central differencing                & \ct{'CENTRAL'}  \\
 Godunov                             & \ct{'GODUNOV'}  \\
 Superbee (VLES, SVLES default)      & \ct{'SUPERBEE'}  \\
-CHARM (DNS, LES default)            & \ct{'CHARM'}  \\ 
+CHARM (DNS, LES default)            & \ct{'CHARM'}  \\
 \hline
 \end{tabular}
 \end{table}
@@ -9961,6 +9960,9 @@ \subsubsection{Limiting the Integration}
 \end{lstlisting}
 would output the total surface area in the volume \ct{XB} where the total heat flux exceeds 10~\unit{kW/m^2}.
 
+\subsubsection{Two-Dimensional and Cylindrical Coordinate Systems}
+
+If the computational domain is two-dimensional or cylindrical, some spatially-integrated quantities are adjusted to eliminate the dependence on the arbitrarily chosen $\delta y$ or $\delta \theta$. For a 2-D domain, a reported volume output will have units of m$^3$/m and an area output will have units of m$^2$/m. For a cylindrical domain, volume and area outputs are reported for the entire cylinder rather than the thin wedge on which the simulation is performed. This might cause confusion in cases where both the input parameters and simulation results involve volumetric or areal quantities; thus, it is good practice to perform a simple test case with a known result to verify that these adjustments have been performed properly.
 
 
 \subsection{Temporally-Integrated Outputs}
@@ -14429,8 +14431,8 @@ \chapter{Error Codes}
 610 \> \ct{HOLE ... Cannot overlap HOLEs with a DEVC or CTRL_ID.} \> Section~\ref{info:HOLE}  \\
 611 \> \ct{OBST ... has a BULK_DENSITY but zero volume.} \> Section~\ref{info:BURN_AWAY}  \\
 612 \> \ct{OBST ... must have a volume to be assigned HT3D.} \> Section~\ref{checkerboard}  \\
-614 \> \ct{OBST_ID ... cannot have a SURF with NODE_ID} \> Section~\ref{info:hvac_geom}  \\    
-615 \> \ct{OBST ... cannot overlap OBST ...}            \> Section~\ref{info:OVERLAY}  \\    
+614 \> \ct{OBST_ID ... cannot have a SURF with NODE_ID} \> Section~\ref{info:hvac_geom}  \\
+615 \> \ct{OBST ... cannot overlap OBST ...}            \> Section~\ref{info:OVERLAY}  \\
     \>                                                  \>                          \\
 701 \> \ct{problem with GEOM, local SURF_ID index ... out of bounds.} \> Section~\ref{info:GEOM_Basics}  \\
 702 \> \ct{problem with GEOM, SURF_IDS not defined properly.} \> Section~\ref{info:first_geom}  \\

Manuals/FDS_Validation_Guide/Burning_Rate_Chapter.tex

@@ -2298,7 +2298,7 @@ \subsection{CSIRO Grassland Fires}
 \subsection{USFS/Catchpole Experiments}
 \label{USFS_Catchpole_Plots}
 
-Figures~\ref{USFS_Catchpole_008} through \ref{USFS_Catchpole_354} present the results of 354 simulations of the USFS/Catchpole experiments. A brief description is given in Sec.~\ref{USFS_Catchpole_Description}. The paper by Catchpole et al.~\cite{Catchpole:CST1998} reports a single rate of spread for each experiment, which is depicted in the figures as a straight black line. The rate of spread of the simulations was calculated by fitting the best line through the data points over a time interval between 10~\% and 90~\% of the observed transit time of the real fire over the 8~m fuel bed. The red dashed line is the best fit line from which the rate of spread is taken.
+Figures~\ref{USFS_Catchpole_008} through \ref{USFS_Catchpole_354} present the results of 354 simulations of the USFS/Catchpole experiments. A brief description is given in Sec.~\ref{USFS_Catchpole_Description}. The paper by Catchpole et al.~\cite{Catchpole:CST1998} reports a single rate of spread for each experiment, which is depicted in the figures as a straight black line. The rate of spread of the simulations was calculated by fitting the best line through the fire position versus time over an interval between 2~m and 7~m from the ignition line, to avoid any edge effects on the 8~m fuel bed. Cases are terminated when the fire extinguishes, so some cases where FDS does not sustain fire spread have few or no observations over this interval.
 
 \newpage
 
@@ -2973,9 +2973,34 @@ \subsection{USFS/Catchpole Experiments}
 \label{USFS_Catchpole_354}
 \end{figure}
 
+\clearpage
+
+With over 300 cases in this study, it can be difficult to evaluate the quality of the FDS predictions as they pertain to the different experimental parameters evaluated. Across the four types of fuel, the input parameters modified were particle surface-to-volume ratio (s), packing ratio (beta), fuel moisture content (M), and wind speed (U). The ratio of the FDS predicted spread rates to the observed spread rates are shown as a function of each of these parameters in Fig.~\ref{USFS_Catchpole_parameter_summary}.
 
+\begin{figure}[ht]
+\begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r}
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_R_v_s} &
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_R_v_beta} \\
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_R_v_M} &
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_R_v_U} \\
+\end{tabular*}
+\caption[USFS/Catchpole, effect of test parameters on spread rate prediction.]{Ratio of FDS to observed spread rate as a function of test matrix parameters for USFS/Catchpole cases. For context, dashed lines represent the region of plus/minus 20~\% error in predicted spread.}
+\label{USFS_Catchpole_parameter_summary}
+\end{figure}
 
+Additionally, neither Fig.~\ref{USFS_Catchpole_parameter_summary} nor the summary scatter plot in Fig.~\ref{RoS_Summary} include information on the cases for which FDS fails to reproduce a spreading fire. All cases tested have a reported spread rate from the experiments, so instances where FDS does not produce sustained fire spread are indicative of limitations in the current model representation. In order to guide future development efforts, Fig.~\ref{USFS_Catchpole_no_spread} attempts to summarize the model performance over the full parameter space. Gray points in the background correspond to successful spread predictions, while colored points in the foreground are those for which FDS does not sustain fire spread. This highlights problematic areas of the test matrix, such as zero wind speed cases with low values of packing ratio.
 
+\begin{figure}[ht]
+\centering
+\begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r}
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_no_spread_MF} &
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_no_spread_PPMC} \\
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_no_spread_EXSC} &
+\includegraphics[height=2.15in]{SCRIPT_FIGURES/USFS_Catchpole/Catchpole_no_spread_EX} \\
+\end{tabular*}
+\caption[USFS/Catchpole, conditions leading to failed spread predictions]{Parameter combinations for all USFS/Catchpole cases. Gray points in the background are cases where FDS successfully predicts spread, and colored points in the foreground represent failed (no spread) predictions. Both the point size and color are scaled according to the packing ratio (beta) value for a given test.}
+\label{USFS_Catchpole_no_spread}
+\end{figure}
 
 \clearpage
 

Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex

@@ -3999,29 +3999,37 @@ \section{Mixing with detailed chemistry (\ct{edc\_mixing\_cvode})}
 
 This implementation has been validated against results obtained using Cantera~\cite{cantera:2023}. Note that Cantera does not provide a direct method to solve the coupled mixing-chemistry ODE system as done here in FDS; therefore, in Cantera an operator splitting approach was employed, where mixing and chemistry are handled sequentially. Specifically, within each CFD timestep (0.1 s), substeps for chemistry calculations in Cantera are kept sufficiently small (0.0001 s) to resolve the fast dynamics. All comparisons are carried out using the GRI-Mech 3.0 chemical mechanism~\cite{gri3:1999}.
 
-Figure~\ref{fig:edc_mixing_cvode_onecfdstep_vary_zeta0} presents comparisons of temperature and OH mass fraction between FDS and Cantera within a single CFD timestep for five constant-volume ignition delay cases, with varying initial unmixed fractions \(\zeta_0 = 1.0, 0.75, 0.5, 0.25, 0.0\).
+Figures~\ref{fig:edc_mixing_cvode_onecfdstep_vary_zeta0_1} and~\ref{fig:edc_mixing_cvode_onecfdstep_vary_zeta0_2} presents comparisons of temperature and OH mass fraction between FDS and Cantera within a single CFD timestep for five constant-volume ignition delay cases, with varying initial unmixed fractions \(\zeta_0 = 1.0, 0.75, 0.5, 0.25, 0.0\).
 
 Figure~\ref{fig:edc_mixing_cvode_multicfdstep_vary_zeta0} shows similar comparisons over multiple CFD timesteps for the same ignition delay cases and initial unmixed fractions. These results also verify that elemental mass is conserved throughout the entire process.
 
 Finally, Figure~\ref{fig:edc_mixing_cvode_onecfdstep_vary_taumix} compares temperature and OH mass fraction for a single CFD timestep across five cases with varying mixing times \(\tau_{\text{mix}} = 0.1, 0.01, 0.001, 0.0001, 0.00001\) s. As the mixing time decreases, the solution becomes independent of \(\tau_{\text{mix}}\), approaching the well-stirred reactor limit corresponding to (\(\zeta_0 = 0.0\)).
 
 \begin{figure}[p]
 \begin{tabular*}{\textwidth}{lr}
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta1p0_TMP} &
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta1p0_OH} \\
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p75_TMP} &
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p75_OH} \\
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p5_TMP} &
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p5_OH} \\
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p25_TMP} &
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p25_OH} \\
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p0_TMP} &
-\includegraphics[height=1.6in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p0_OH} 
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta1p0_TMP} &
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta1p0_OH} \\
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p75_TMP} &
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p75_OH} \\
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p5_TMP} &
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p5_OH} \\
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p25_TMP} &
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p25_OH}
 \end{tabular*}
 \caption[Results of the \ct{edc\_mixing\_cvode} test cases]{Comparison of CVODE substeps for a single CFD step of 0.1 s, showing the effect of varying the initial unmixed fraction between FDS and Cantera. The simulation uses the Methane GRI mechanism with an equivalence ratio of 0.6, a mixing time of 0.01 s, and an initial temperature of 1200 K.}
-\label{fig:edc_mixing_cvode_onecfdstep_vary_zeta0}
+\label{fig:edc_mixing_cvode_onecfdstep_vary_zeta0_1}
 \end{figure}
 
+\begin{figure}[p]
+\begin{tabular*}{\textwidth}{lr}
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p0_TMP} &
+\includegraphics[height=2.2in]{SCRIPT_FIGURES/EDC_OneCFDStep_Methane_grimech30_Zeta0p0_OH}
+\end{tabular*}
+\caption[Results of the \ct{edc\_mixing\_cvode} test cases]{Continued from Figure~\ref{fig:edc_mixing_cvode_onecfdstep_vary_zeta0_1}.}
+\label{fig:edc_mixing_cvode_onecfdstep_vary_zeta0_2}
+\end{figure}
+
+
 \begin{figure}[p]
 \begin{tabular*}{\textwidth}{lcr}
 \includegraphics[height=1.6in,width=0.32\textwidth]{SCRIPT_FIGURES/EDC_MultiCFDStep_Methane_grimech30_Zeta1p0_TMP} &
@@ -5141,7 +5149,7 @@ \subsection{Case 8: 2-D Heat Transfer with Temperature-Dependent Conductivity}
 \FloatBarrier
 
 
-\subsection{Case 9: 2-D Heat Transfer in a Composite Section with Temperature-Dependent Conductivity}
+\subsection{Case 9: 2-D Heat Transfer in a Composite Section, Variable Conductivity}
 \label{SFPE_Case_9}
 
 A hollow square metal tube ($\rho=7850$~kg/m$^3$, $c=0.6$~kJ/(kg~K), $\epsilon=0.8$) is filled with an insulation material ($k=0.05$~W/(m~K), $\rho=50$~kg/m$^3$, $c=1$~kJ/(kg~K)). The thermal conductivity of the metal tube varies linearly with temperature such that its value is 54.7~W/(m~K) at 0~°C, 27.3~W/(m~K) at 800~°C, and 27.3~W/(m~K) at 1200~°C. The tube walls are 0.5~mm thick, and the exterior width of the assembly is 0.201~m. The surrounding air temperature is 1000~°C, and the initial temperature of the assembly is 0~°C. Assuming that the heating is by convection and radiation, and that the heat transfer coefficient is 10~W/(m$^2$~K), calculate the temperature at the center of the tube as a function of time (Fig.~\ref{fig:SFPE_Case_9}).

Source/cons.f90

@@ -295,6 +295,7 @@ MODULE GLOBAL_CONSTANTS
 LOGICAL :: STORE_FIRE_RESIDENCE=.FALSE.             !< Flag for tracking residence time of spreading fire front over a surface
 LOGICAL :: STORE_LS_SPREAD_RATE=.FALSE.             !< Flag for outputting local level set spread rate magnitude
 LOGICAL :: TEST_NEW_CHAR_MODEL=.FALSE.              !< Flag to envoke new char model
+LOGICAL :: FLUX_LIMITER_MW_CORRECTION=.TRUE.        !< Flag for MW correction ensure consistent equation of state at face
 
 INTEGER, ALLOCATABLE, DIMENSION(:) :: CHANGE_TIME_STEP_INDEX      !< Flag to indicate if a mesh needs to change time step
 INTEGER, ALLOCATABLE, DIMENSION(:) :: SETUP_PRESSURE_ZONES_INDEX  !< Flag to indicate if a mesh needs to keep searching for ZONEs

Source/data.f90

@@ -1360,42 +1360,6 @@ SUBROUTINE DEFINE_OUTPUT_QUANTITIES
 OUTPUT_QUANTITY(553)%UNITS  = 'm/s'
 OUTPUT_QUANTITY(553)%SHORT_NAME  = 'V_LS'
 
-OUTPUT_QUANTITY(560)%NAME  = 'BFX'
-OUTPUT_QUANTITY(560)%UNITS  = ''
-OUTPUT_QUANTITY(560)%SHORT_NAME  = 'bfx'
-OUTPUT_QUANTITY(560)%CELL_POSITION = CELL_FACE
-OUTPUT_QUANTITY(560)%IOR = 1
-
-OUTPUT_QUANTITY(561)%NAME  = 'BFY'
-OUTPUT_QUANTITY(561)%UNITS  = ''
-OUTPUT_QUANTITY(561)%SHORT_NAME  = 'bfy'
-OUTPUT_QUANTITY(561)%CELL_POSITION = CELL_FACE
-OUTPUT_QUANTITY(561)%IOR = 2
-
-OUTPUT_QUANTITY(562)%NAME  = 'BFZ'
-OUTPUT_QUANTITY(562)%UNITS  = ''
-OUTPUT_QUANTITY(562)%SHORT_NAME  = 'bfz'
-OUTPUT_QUANTITY(562)%CELL_POSITION = CELL_FACE
-OUTPUT_QUANTITY(562)%IOR = 3
-
-OUTPUT_QUANTITY(563)%NAME  = 'BFX MINUS'
-OUTPUT_QUANTITY(563)%UNITS  = ''
-OUTPUT_QUANTITY(563)%SHORT_NAME  = 'bfx-'
-OUTPUT_QUANTITY(563)%CELL_POSITION = CELL_FACE
-OUTPUT_QUANTITY(563)%IOR = 1
-
-OUTPUT_QUANTITY(564)%NAME  = 'BFY MINUS'
-OUTPUT_QUANTITY(564)%UNITS  = ''
-OUTPUT_QUANTITY(564)%SHORT_NAME  = 'bfy-'
-OUTPUT_QUANTITY(564)%CELL_POSITION = CELL_FACE
-OUTPUT_QUANTITY(564)%IOR = 2
-
-OUTPUT_QUANTITY(565)%NAME  = 'BFZ MINUS'
-OUTPUT_QUANTITY(565)%UNITS  = ''
-OUTPUT_QUANTITY(565)%SHORT_NAME  = 'bfz-'
-OUTPUT_QUANTITY(565)%CELL_POSITION = CELL_FACE
-OUTPUT_QUANTITY(565)%IOR = 3
-
 ! Boundary Quantities (Negative indices)
 
 OUTPUT_QUANTITY(-1)%NAME = 'RADIATIVE HEAT FLUX'