diff --git a/Manuals/Copy_Firebot_FiguresGH.sh b/Manuals/Copy_Firebot_FiguresGH.sh index ff89ae9caac..bf2be4772c6 100755 --- a/Manuals/Copy_Firebot_FiguresGH.sh +++ b/Manuals/Copy_Firebot_FiguresGH.sh @@ -1,7 +1,12 @@ #!/bin/bash +LATEST=$1 +if [ "$LATEST" != "" ]; then + LATEST=_latest +fi + DOWNLOADFIGURES () { - FILE=$1_figures + FILE=$1_figures$LATEST echo ***Downloading $FILE.tar.gz gh release download FDS_TEST -p $FILE.tar.gz -R github.com/firemodels/test_bundles --clobber if [ ! -e $FILE.tar.gz ]; then @@ -31,25 +36,25 @@ FBVG=$BASEDIR/FIGS/FDS_VERG FBVAL=$BASEDIR/FIGS/FDS_VALG # Copy Tech Guide Figures -if [ -e $BASEDIR/FDS_TG_figures.tar.gz ]; then +if [ -e $BASEDIR/FDS_TG_figures$LATEST.tar.gz ]; then cd $FBTG - tar xf $BASEDIR/FDS_TG_figures.tar.gz + tar xf $BASEDIR/FDS_TG_figures$LATEST.tar.gz cp * $BASEDIR/FDS_Technical_Reference_Guide/SCRIPT_FIGURES/ echo FDS Technical Guide figures copied to FDS_Technical_Reference_Guide/SCRIPT_FIGURES fi # Copy User's Guide Figures -if [ -e $BASEDIR/FDS_UG_figures.tar.gz ]; then +if [ -e $BASEDIR/FDS_UG_figures$LATEST.tar.gz ]; then cd $FBUG - tar xf $BASEDIR/FDS_UG_figures.tar.gz + tar xf $BASEDIR/FDS_UG_figures$LATEST.tar.gz cp * $BASEDIR/FDS_User_Guide/SCRIPT_FIGURES/ echo FDS Users Guide figures copied to FDS_User_Reference_Guide/SCRIPT_FIGURES fi # Copy Verification Guide Figures -if [ -e $BASEDIR/FDS_VERG_figures.tar.gz ]; then +if [ -e $BASEDIR/FDS_VERG_figures$LATEST.tar.gz ]; then cd $FBVG - tar xf $BASEDIR/FDS_VERG_figures.tar.gz + tar xf $BASEDIR/FDS_VERG_figures$LATEST.tar.gz cp *.pdf $BASEDIR/FDS_Verification_Guide/SCRIPT_FIGURES/. cp *.png $BASEDIR/FDS_Verification_Guide/SCRIPT_FIGURES/. cp *.tex $BASEDIR/FDS_Verification_Guide/SCRIPT_FIGURES/. @@ -58,9 +63,9 @@ if [ -e $BASEDIR/FDS_VERG_figures.tar.gz ]; then fi # Copy Validation Guide Figures -if [ -e $BASEDIR/FDS_VALG_figures.tar.gz ]; then +if [ -e $BASEDIR/FDS_VALG_figures$LATEST.tar.gz ]; then cd $FBVAL - tar xf $BASEDIR/FDS_VALG_figures.tar.gz + tar xf $BASEDIR/FDS_VALG_figures$LATEST.tar.gz cp -R * $BASEDIR/FDS_Validation_Guide/SCRIPT_FIGURES echo FDS Validation Guide Figures copied to FDS_Validation_Guide_Reference_Guide/SCRIPT_FIGURES fi diff --git a/Manuals/FDS_User_Guide/FDS_User_Guide.tex b/Manuals/FDS_User_Guide/FDS_User_Guide.tex index 66c6a2d2439..e875ed3ea64 100644 --- a/Manuals/FDS_User_Guide/FDS_User_Guide.tex +++ b/Manuals/FDS_User_Guide/FDS_User_Guide.tex @@ -4621,7 +4621,7 @@ \subsubsection{Enthalpy} \be c_{p,\alpha} = \frac{\gamma}{\gamma-1} \frac{\R}{W_\alpha} \ee -The ratio of specific heats, \ct{GAMMA}, is 1.4 by default and can be changed on the \ct{MISC} line. If you want all the gas specific heats to follow this relation, set \ct{CONSTANT_SPECIFIC_HEAT_RATIO=T} on the \ct{MISC} line. For high molecular weight species, use of the default gamma will result in very low values of the specific heat which can cause issues with the default extinction model. In this case it is recommended that the specific heat be defined. If the FDS is determining the specific heat using this relation, then it is recommended that you check the \ct{CHID.out} and verify that reasonable specific heat values have been created. +The ratio of specific heats, \ct{GAMMA}, is 1.4 by default and can be changed on the \ct{MISC} line. If you want all the gas specific heats to follow this relation, set \ct{CONSTANT_SPECIFIC_HEAT_RATIO=T} on the \ct{MISC} line. For high molecular weight species, use of the default gamma will result in very low values of the specific heat. This can cause issues with temperature in regions with high fuel mass fractions and can cause issues with the \ct{EXTINCTION 2} extinction model. For high molecular weight species it is recommended that the specific heat be defined. If the FDS is determining the specific heat using \ct{GAMMA}, then it is recommended that you check the \ct{CHID.out} file and verify that reasonable specific heat values have been created. Note that species used in chemical reactions must either be predefined or have an explicitly defined $h(T_{\rm ref})$. The exception is the species defined as \ct{FUEL} on the \ct{REAC} input as long as a \ct{HEAT_OF_COMBUSTION} is defined or the reaction is a simple chemistry reaction where \ct{EPUMO2} applies. Specifying either \ct{SPECIFIC_HEAT} or \ct{RAMP_CP} is considered to have explicitly defined $h(T_{\rm ref})$. If you still wish to use the default \ct{EPUMO2} or 13,100~kJ/kg, then you will need to explicitly define it on the \ct{REAC} input. @@ -6888,12 +6888,12 @@ \subsection{Basic Equations} \begin{figure}[p] \begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r} - \includegraphics[width=3.2in]{FIGURES/vel_L=-100} & - \includegraphics[width=3.2in]{FIGURES/tmp_L=-100} \\ - \includegraphics[width=3.2in]{FIGURES/vel_L=100} & - \includegraphics[width=3.2in]{FIGURES/tmp_L=100} \\ - \includegraphics[width=3.2in]{FIGURES/vel_L=-350} & - \includegraphics[width=3.2in]{FIGURES/tmp_L=-350} + \includegraphics[width=3.2in]{FIGURES/Monin_Obukhov_vel_L_-100} & + \includegraphics[width=3.2in]{FIGURES/Monin_Obukhov_tmp_L_-100} \\ + \includegraphics[width=3.2in]{FIGURES/Monin_Obukhov_vel_L_100} & + \includegraphics[width=3.2in]{FIGURES/Monin_Obukhov_tmp_L_100} \\ + \includegraphics[width=3.2in]{FIGURES/Monin_Obukhov_vel_L_-350} & + \includegraphics[width=3.2in]{FIGURES/Monin_Obukhov_tmp_L_-350} \end{tabular*} \caption[Sample vertical wind and temperature profiles]{Sample vertical wind and temperature profiles.} \label{wind_profiles} @@ -12267,7 +12267,7 @@ \section{\texorpdfstring{{\tt DEVC}}{DEVC} (Device Parameters)} \ct{SMOOTHING_TIME} & Real & Section~\ref{info:basic_control} & s & \\ \hline \ct{SPATIAL_STATISTIC} & Character & Section~\ref{info:statistics} & & \\ \hline \ct{SPEC_ID} & Character & Section~\ref{info:gasoutputquantities} & & \\ \hline -\ct{STATISTICS_END} & Real & Section~\ref{info:rmscovcorr} & s & \ct{T_BEGIN}\\ \hline +\ct{STATISTICS_END} & Real & Section~\ref{info:rmscovcorr} & s & \ct{T_END}\\ \hline \ct{STATISTICS_START} & Real & Section~\ref{info:rmscovcorr} & s & \ct{T_BEGIN}\\ \hline \ct{SURF_ID} & Character & Section~\ref{info:statistics} & & \\ \hline \ct{TEMPORAL_STATISTIC} & Character & Section~\ref{info:statistics} & & \\ \hline diff --git a/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_Profiles.py b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_Profiles.py new file mode 100644 index 00000000000..ca28fcd3389 --- /dev/null +++ b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_Profiles.py @@ -0,0 +1,83 @@ + +# Atmospheric Boundary Layer profiles, based on M-O theory as described in +# "Falcon Series Data Report", GRI-89/0138, June 1990. +# These plots are used as illustrations in the FDS User's Guide. + +import numpy as np +import matplotlib.pyplot as plt +import fdsplotlib + +plot_style = fdsplotlib.get_plot_style('fds') + +rho_0 = 1.2 +g = 9.81 +cp = 1005 +kappa = 0.41 +p_0 = 100000 + +L = np.array([-100,100,-350]) +Ls = np.array(['-100','100','-350']) +z_0 = np.array([0.01,0.005,0.005]) +z_r = 2.0 +u_r = 5.0 + +for i in range(3): + + p_r = p_0 - rho_0 * g * (z_r - z_0[i]) + T_r = 20 + 273.15 + theta_r = T_r * (p_0 / p_r) ** 0.285 + + if L[i] < 0: + x_r = (1 - 16 * z_r / L[i]) ** 0.25 + psi_m_r = 2 * np.log((1 + x_r) / 2) + np.log((1 + x_r**2) / 2) - 2 * np.arctan(x_r) + np.pi / 2 + psi_h_r = 2 * np.log((1 + x_r**2) / 2) + else: + psi_m_r = -5 * z_r / L[i] + psi_h_r = psi_m_r + + u_star = kappa * u_r / (np.log(z_r / z_0[i]) - psi_m_r) + theta_0 = theta_r / (1 + u_star**2 * (np.log(z_r / z_0[i]) - psi_h_r) / (g * kappa**2 * L[i])) + theta_star = u_star**2 * theta_0 / (g * kappa * L[i]) + + z = np.array([z_0[i], 10*z_0[i], 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 25, 30, 50, 100]) + + # Figure 1 - Velocity profile + + if L[i] < 0: + x = (1 - 16 * z / L[i]) ** 0.25 + psi_m = 2 * np.log((1 + x) / 2) + np.log((1 + x**2) / 2) - 2 * np.arctan(x) + np.pi / 2 + psi_h = 2 * np.log((1 + x**2) / 2) + else: + psi_m = -5 * z / L[i] + psi_h = psi_m + + u = (u_star / kappa) * (np.log(z / z_0[i]) - psi_m) + theta = theta_0 + (theta_star / kappa) * (np.log(z / z_0[i]) - psi_h) + T = theta * (p_0 / (p_0 - rho_0 * g * (z - z_0[i]))) ** -0.285 + + fig = fdsplotlib.plot_to_fig(x_data=u, y_data=z, marker_style='ko-', + x_label='Velocity (m/s)', y_label='Height (m)', + x_min=0, x_max=10, y_min=0, y_max=30) + + ax1 = fig.axes[0] + ax1.text(0.05, 0.90, r'$u(2 \; \mathrm{m})=5$ m/s', transform=ax1.transAxes, verticalalignment='top') + ax1.text(0.05, 0.80, f'$L={L[i]:4.0f}$ m', transform=ax1.transAxes, verticalalignment='top') + ax1.text(0.05, 0.70, f'$z_0={z_0[i]:5.3f}$ m', transform=ax1.transAxes, verticalalignment='top') + + plt.savefig(f'Monin_Obukhov_vel_L_{Ls[i]}.pdf', format='pdf') + plt.close() + + # Figure 2 - Temperature profile + + fig = fdsplotlib.plot_to_fig(x_data=T-273.15, y_data=z, marker_style='ko-', + x_label='Temperature (°C)', y_label='Height (m)', + x_min=15, x_max=25, y_min=0, y_max=30) + + ax2 = fig.axes[0] + ax2.text(0.05, 0.90, r'$T(2 \; \mathrm{m})=20\;^\circ$C', transform=ax2.transAxes, verticalalignment='top') + ax2.text(0.05, 0.80, f'$L={L[i]:4.0f}$ m', transform=ax2.transAxes, verticalalignment='top') + ax2.text(0.05, 0.70, f'$z_0={z_0[i]:5.3f}$ m', transform=ax2.transAxes, verticalalignment='top') + + plt.savefig(f'Monin_Obukhov_tmp_L_{Ls[i]}.pdf', format='pdf') + plt.close() + diff --git a/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_-100.pdf b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_-100.pdf new file mode 100644 index 00000000000..cdf5e2e33b9 Binary files /dev/null and b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_-100.pdf differ diff --git a/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_-350.pdf b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_-350.pdf new file mode 100644 index 00000000000..cc30142648f Binary files /dev/null and b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_-350.pdf differ diff --git a/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_100.pdf b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_100.pdf new file mode 100644 index 00000000000..40ab7937aa5 Binary files /dev/null and b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_tmp_L_100.pdf differ diff --git a/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_-100.pdf b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_-100.pdf new file mode 100644 index 00000000000..5b548fce95a Binary files /dev/null and b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_-100.pdf differ diff --git a/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_-350.pdf b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_-350.pdf new file mode 100644 index 00000000000..e53857b0462 Binary files /dev/null and b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_-350.pdf differ diff --git a/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_100.pdf b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_100.pdf new file mode 100644 index 00000000000..b59935fdfae Binary files /dev/null and b/Manuals/FDS_User_Guide/FIGURES/Monin_Obukhov_vel_L_100.pdf differ diff --git a/Utilities/Python/scripts/spray_pattern.py b/Manuals/FDS_User_Guide/FIGURES/spray_pattern.py similarity index 88% rename from Utilities/Python/scripts/spray_pattern.py rename to Manuals/FDS_User_Guide/FIGURES/spray_pattern.py index b97fef75dab..2a85948e4c2 100644 --- a/Utilities/Python/scripts/spray_pattern.py +++ b/Manuals/FDS_User_Guide/FIGURES/spray_pattern.py @@ -7,8 +7,6 @@ plot_style = fdsplotlib.get_plot_style('fds') -pltdir = '../../Manuals/FDS_User_Guide/FIGURES/' - b = [0, 5, 10, 100, 1000] mu = np.array([0, 10, 20]) * np.pi / 180 mu_str = ['0', '10', '20'] @@ -34,6 +32,6 @@ ax.invert_yaxis() ax.text(45, 0.1, r'$\mu=$' + str(mu[j]*180/np.pi) + r', $\phi_{\rm max}=$' + str(phi_max[k]*180/np.pi) + ' degrees', fontsize=plot_style['Key_Font_Size']) - plt.savefig(pltdir + 'spray_pattern_mu_' + mu_str[j] + '_phimax_' + phi_str[k] + '.pdf', format='pdf') + plt.savefig('spray_pattern_mu_' + mu_str[j] + '_phimax_' + phi_str[k] + '.pdf', format='pdf') plt.close() diff --git a/Manuals/FDS_User_Guide/FIGURES/tmp_L=-100.pdf b/Manuals/FDS_User_Guide/FIGURES/tmp_L=-100.pdf deleted file mode 100644 index 2a970edf19e..00000000000 Binary files a/Manuals/FDS_User_Guide/FIGURES/tmp_L=-100.pdf and /dev/null differ diff --git a/Manuals/FDS_User_Guide/FIGURES/tmp_L=-350.pdf b/Manuals/FDS_User_Guide/FIGURES/tmp_L=-350.pdf deleted file mode 100644 index bd37f6bc6e7..00000000000 Binary files a/Manuals/FDS_User_Guide/FIGURES/tmp_L=-350.pdf and /dev/null differ diff --git a/Manuals/FDS_User_Guide/FIGURES/tmp_L=100.pdf b/Manuals/FDS_User_Guide/FIGURES/tmp_L=100.pdf deleted file mode 100644 index 204f69646c8..00000000000 Binary files a/Manuals/FDS_User_Guide/FIGURES/tmp_L=100.pdf and /dev/null differ diff --git a/Manuals/FDS_User_Guide/FIGURES/vel_L=-100.pdf b/Manuals/FDS_User_Guide/FIGURES/vel_L=-100.pdf deleted file mode 100644 index b16733d9934..00000000000 Binary files a/Manuals/FDS_User_Guide/FIGURES/vel_L=-100.pdf and /dev/null differ diff --git a/Manuals/FDS_User_Guide/FIGURES/vel_L=-350.pdf b/Manuals/FDS_User_Guide/FIGURES/vel_L=-350.pdf deleted file mode 100644 index 03c5a6ffc4a..00000000000 Binary files a/Manuals/FDS_User_Guide/FIGURES/vel_L=-350.pdf and /dev/null differ diff --git a/Manuals/FDS_User_Guide/FIGURES/vel_L=100.pdf b/Manuals/FDS_User_Guide/FIGURES/vel_L=100.pdf deleted file mode 100644 index bfa7c97bcde..00000000000 Binary files a/Manuals/FDS_User_Guide/FIGURES/vel_L=100.pdf and /dev/null differ diff --git a/Source/part.f90 b/Source/part.f90 index b2e0872b197..1dde6c46582 100644 --- a/Source/part.f90 +++ b/Source/part.f90 @@ -1013,9 +1013,9 @@ SUBROUTINE INSERT_VOLUMETRIC_PARTICLES ! If the INIT volume is outside the current mesh, return IF (IN%SHAPE/='RING' .AND. IN%SHAPE/='LINE') THEN - IF ((IN_X1-XF)>50._EB*TWO_EPSILON_EB .OR. (IN_X2-XS)<-50._EB*TWO_EPSILON_EB .OR. & - (IN_Y1-YF)>50._EB*TWO_EPSILON_EB .OR. (IN_Y2-YS)<-50._EB*TWO_EPSILON_EB .OR. & - (IN_Z1-ZF)>50._EB*TWO_EPSILON_EB .OR. (IN_Z2-ZS)<-50._EB*TWO_EPSILON_EB) RETURN + IF ((IN_X1-XF)>-50._EB*TWO_EPSILON_EB .OR. (IN_X2-XS)<50._EB*TWO_EPSILON_EB .OR. & + (IN_Y1-YF)>-50._EB*TWO_EPSILON_EB .OR. (IN_Y2-YS)<50._EB*TWO_EPSILON_EB .OR. & + (IN_Z1-ZF)>-50._EB*TWO_EPSILON_EB .OR. (IN_Z2-ZS)<50._EB*TWO_EPSILON_EB) RETURN ELSEIF (IN%SHAPE=='RING') THEN IF (RING_MESH_INTERSECTION_ARC(NM,IN%X0,IN%Y0,IN%RADIUS) 1) CALL MPI_ALLREDUCE(TRN_ME(1),TRN_ME(2),1,MPI_INTEGER,MPI_SUM,MPI_COMM_WORLD,IERR) IF (TRN_ME(2) > 0) THEN ! Meshes at different refinement levels. Not Unsupported. - IF (MY_RANK == 0) WRITE(LU_ERR,*) 'GLMAT Setup Error : Meshes at different refinement levels unsupported.' + IF (MY_RANK == 0) WRITE(LU_ERR,*) 'GLMAT Setup Error: Meshes at different refinement levels currently unsupported for &GEOM.' SUPPORTED_MESH = .FALSE. STOP_STATUS = SETUP_STOP RETURN diff --git a/Source/read.f90 b/Source/read.f90 index ff78279a215..382476f5c76 100644 --- a/Source/read.f90 +++ b/Source/read.f90 @@ -12999,12 +12999,12 @@ SUBROUTINE READ_INIT MESH_LOOP: DO NM=1,NMESHES M=>MESHES(NM) IF (XYZ(1)>=M%XS .AND. XYZ(1)<=M%XF .AND. XYZ(2)>=M%YS .AND. XYZ(2)<=M%YF .AND. XYZ(3)>=M%ZS .AND. XYZ(3)<=M%ZF) THEN - IF (ABS(XYZ(1)-M%XS)tol - %--------- - %NOTE: coefficients were found from NIST Webbook - %--------- - - %initial temperature coeffs - coeff0(:,1) = [28.98641;1.853978;-9.647459;16.63537;0.000117;-8.671914;0.0]; - coeff0(:,2) = [-0.703029;108.4773;-42.52157;5.862788;0.678565;-76.84376;-74.87310]; - coeff0(:,3) = [31.32234;-20.23531;57.86644;-36.50624;-0.007374;-8.903471;0.0]; - coeff0(:,4) = [24.99735;55.18696;-33.69137;7.948387;-0.136638;-403.6075;-393.5224]; - coeff0(:,5) = [30.09200;6.832514;6.793435;-2.534480;0.082139;-250.8810;-241.8264]; - - %nitrogen cp coeffs [J/mol K] - if Tf_guess <=500 - coeff(:,1) = [28.98641;1.853978;-9.647459;16.63537;0.000117;-8.671914;0.0]; - elseif 500 < Tf_guess <= 2000 - coeff(:,1) = [19.50583;19.88705;-8.598535;1.369784;0.527601;-4.935202;0.0]; - else - coeff(:,1) = [35.51872;1.128728;-0.196103;0.014662;-4.553760;-18.97091;0.0]; - end - %methane cp coeffs [J/mol K] - if Tf_guess <= 1300 - coeff(:,2) = [-0.703029;108.4773;-42.52157;5.862788;0.678565;-76.84376;-74.87310]; - else - coeff(:,2) = [85.81217;11.26467;-2.114146;0.138190;-26.42221;-153.5327;-74.87310]; - end - %oxygen cp coeffs [J/mol K] - if Tf_guess <=700 - coeff(:,3) = [31.32234;-20.23531;57.86644;-36.50624;-0.007374;-8.903471;0.0]; - elseif 700 < Tf_guess <= 2000 - coeff(:,3) = [30.03235;8.772972;-3.988133;0.788313;-0.741599;-11.32468;0.0]; - else - coeff(:,3) = [20.91111;10.72071;-2.020498;0.146449;9.245722;5.337651;0.0]; - end - %carbon dooxide cp coeffs [J/mol K] - if Tf_guess <= 1200 - coeff(:,4) = [24.99735;55.18696;-33.69137;7.948387;-0.136638;-403.6075;-393.5224]; - else - coeff(:,4) = [58.16639;2.720074;-0.492289;0.038844;-6.447293;-425.9186;-393.5224]; - end - %water vapor cp coeffs [J/mol K] - if Tf_guess <= 1700 - coeff(:,5) = [30.09200;6.832514;6.793435;-2.534480;0.082139;-250.8810;-241.8264]; - else - coeff(:,5) = [41.96246;8.622053;-1.499780;0.098199;-11.15764;-272.1797;-241.8264]; - end - - t=Tf_guess./1000; - t0=temp_0/1000; - for j=1:3 - for i=1:5 - del_h(i,j) = ((coeff(1,i)*t(j) + (1/2)*coeff(2,i)*t(j)^2 + (1/3)*coeff(3,i)*t(j)^3 + (1/4)*coeff(4,i)*t(j)^4 - coeff(5,i)/t(j) + coeff(6,i) - coeff(7,i))... - -(coeff0(1,i)*t0 + (1/2)*coeff0(2,i)*t0^2 + (1/3)*coeff0(3,i)*t0^3 + (1/4)*coeff0(4,i)*t0^4 - coeff0(5,i)/t0 + coeff0(6,i) - coeff0(7,i)))... - * 1000; - end - end - - h_fc = Nf(4)*y_hf(4)+Nf(5)*y_hf(5)-N0(2)*y_hf(2); - del_hN = Nf(1).*del_h(1,:)+Nf(4).*del_h(4,:)+Nf(5).*del_h(5,:); - RT0 = sum(N0)*R*temp_0; - RTf = -sum(Nf)*R*Tf_guess; - - Tf_solve = h_fc + del_hN + RT0 + RTf'; - tol0 = Tf_guess(1) - Tf_guess(3); - - if Tf_solve(2) < 0 - if Tf_solve(1) > 0 - mp = 0.5*(Tf_guess(1)+Tf_guess(2)); - Tf_guess = [Tf_guess(1);mp;Tf_guess(2)]; - else - mp = 0.5*(Tf_guess(2)+Tf_guess(3)); - Tf_guess = [Tf_guess(2);mp;Tf_guess(3)]; - end - else - if Tf_solve(1) < 0 - mp = 0.5*(Tf_guess(1)+Tf_guess(2)); - Tf_guess = [Tf_guess(1);mp;Tf_guess(2)]; - else - mp = 0.5*(Tf_guess(2)+Tf_guess(3)); - Tf_guess = [Tf_guess(2);mp;Tf_guess(3)]; - end - end - - count = count+1; -end -Tf=Tf_guess(2); -Pf=(Tf_guess(2)*sum(Nf)*R)/volume; -dP = Pf-pres_0; -Tf = Tf - 273.15; - -%----------------------- -% Create SS Vector of T and P -%----------------------- - -tss=[1;1.25;1.5;1.75;2]; -Tf=[Tf;Tf;Tf;Tf;Tf]; -dP=[dP;dP;dP;dP;dP]; - -%------------------------------------------------ -% Write Expected Data CSV File -%------------------------------------------------ -yff(:,1) = tspan; - -header1 = {'Time','O2','CH4','CO2','H2O'}; -filename1 = '../../Verification/Species/reactionrate_EDC_flim_1step_CH4_spec.csv'; -fid = fopen(filename1,'wt'); -fprintf(fid,'%s, %s, %s, %s, %s\n',header1{:}); -for j=1:length(tspan) - fprintf(fid,'%f, %f, %f, %f, %f\n',yff(j,:)); -end -fclose(fid); - -header1 = {'Time','TEMP','PRES'}; -filename1 = '../../Verification/Species/reactionrate_EDC_flim_1step_CH4_temppres.csv'; -fid = fopen(filename1,'wt'); -fprintf(fid,'%s, %s, %s\n',header1{:}); -for j=1:length(tss) - fprintf(fid,'%f, %f, %f\n',tss(j),Tf(j),dP(j)); -end -fclose(fid); diff --git a/Utilities/Matlab/scripts/EDC_methane_nonmix.m b/Utilities/Matlab/scripts/EDC_methane_nonmix.m deleted file mode 100644 index f47dfa92010..00000000000 --- a/Utilities/Matlab/scripts/EDC_methane_nonmix.m +++ /dev/null @@ -1,199 +0,0 @@ -%----------------------- -% C Weinschenk -% Verification of EDM (Species,Temp,Pres) -% Fuel: Methane -% Non premix case/fuel rich -% 1/2012 -%----------------------- - -clear all -close all - -%----------------------- -% Initialize species -%----------------------- - -%------------- -% vector key -% (1) = nitrogen -% (2) = methane -% (3) = oxygen -% (4) = carbon dioxide -% (5) = water vapor -%------------- - -temp_0 = 298.15; % intial temperature [K] -volume = 0.001; % volume [m^3] -rho=1.1294838; % density [kg/m^3] -R = 8.3145; % gas constnat [J/mol K] -mass = rho*volume; % mass [kg] - -yf0 = [0.7252; 0.0582; 0.2166; 0.0; 0.0]; -y_MW = [28.0134; 16.042460; 31.9988; 44.0095; 18.01528]; % [g/mol] - -y_hf = [0.0; -74873; 0.0; -393522; -241826]; % [J/mol] - -N0 = (1000*volume*rho*yf0)./(y_MW); % initial moles - -pres_0 = temp_0*sum(N0)*R/volume; % initial pressure - -%----------------------- -% Setup time vector for integration -%----------------------- -tspan=0:5:60; - -%----------------------- -% Setup vector of moles for integrator -%----------------------- -y0=[N0(1) N0(2) N0(3) N0(4) N0(5)]; - -%----------------------- -% Pass information to integrator -%----------------------- -options = odeset('RelTol',1e-10,'AbsTol',1e-10); -[T,Y]=ode113(@EDC_spec_methane,tspan,y0, options); - -Nf = Y(end,:); - -%----------------------- -% Convert back to mass fraction -%----------------------- -yff(:,1)=(y_MW(1)*Y(:,1))/(1000*rho*volume); -yff(:,2)=(y_MW(2)*Y(:,2))/(1000*rho*volume); -yff(:,3)=(y_MW(3)*Y(:,3))/(1000*rho*volume); -yff(:,4)=(y_MW(4)*Y(:,4))/(1000*rho*volume); -yff(:,5)=(y_MW(5)*Y(:,5))/(1000*rho*volume); - -%----------------------- -% Determine final temperature and pressure -%----------------------- -tol0 = 1e5; -tol = 0.1; -Tf_guess = [1000;2000;3000]; % final temperature guess [K] -count = 0; -while abs(tol0)>tol - %--------- - %NOTE: coefficients were found from NIST Webbook - %--------- - - %initial temperature coeffs - coeff0(:,1) = [28.98641;1.853978;-9.647459;16.63537;0.000117;-8.671914;0.0]; - coeff0(:,2) = [-0.703029;108.4773;-42.52157;5.862788;0.678565;-76.84376;-74.87310]; - coeff0(:,3) = [31.32234;-20.23531;57.86644;-36.50624;-0.007374;-8.903471;0.0]; - coeff0(:,4) = [24.99735;55.18696;-33.69137;7.948387;-0.136638;-403.6075;-393.5224]; - coeff0(:,5) = [30.09200;6.832514;6.793435;-2.534480;0.082139;-250.8810;-241.8264]; - - %nitrogen cp coeffs [J/mol K] - if Tf_guess <=500 - coeff(:,1) = [28.98641;1.853978;-9.647459;16.63537;0.000117;-8.671914;0.0]; - elseif 500 < Tf_guess <= 2000 - coeff(:,1) = [19.50583;19.88705;-8.598535;1.369784;0.527601;-4.935202;0.0]; - else - coeff(:,1) = [35.51872;1.128728;-0.196103;0.014662;-4.553760;-18.97091;0.0]; - end - %methane cp coeffs [J/mol K] - if Tf_guess <= 1300 - coeff(:,2) = [-0.703029;108.4773;-42.52157;5.862788;0.678565;-76.84376;-74.87310]; - else - coeff(:,2) = [85.81217;11.26467;-2.114146;0.138190;-26.42221;-153.5327;-74.87310]; - end - %oxygen cp coeffs [J/mol K] - if Tf_guess <=700 - coeff(:,3) = [31.32234;-20.23531;57.86644;-36.50624;-0.007374;-8.903471;0.0]; - elseif 700 < Tf_guess <= 2000 - coeff(:,3) = [30.03235;8.772972;-3.988133;0.788313;-0.741599;-11.32468;0.0]; - else - coeff(:,3) = [20.91111;10.72071;-2.020498;0.146449;9.245722;5.337651;0.0]; - end - %carbon dioxide cp coeffs [J/mol K] - if Tf_guess <= 1200 - coeff(:,4) = [24.99735;55.18696;-33.69137;7.948387;-0.136638;-403.6075;-393.5224]; - else - coeff(:,4) = [58.16639;2.720074;-0.492289;0.038844;-6.447293;-425.9186;-393.5224]; - end - %water vapor cp coeffs [J/mol K] - if Tf_guess <= 1700 - coeff(:,5) = [30.09200;6.832514;6.793435;-2.534480;0.082139;-250.8810;-241.8264]; - else - coeff(:,5) = [41.96246;8.622053;-1.499780;0.098199;-11.15764;-272.1797;-241.8264]; - end - - t=Tf_guess./1000; - t0=temp_0/1000; - for j=1:3 - for i=1:5 - del_h(i,j) = ((coeff(1,i)*t(j) + (1/2)*coeff(2,i)*t(j)^2 + (1/3)*coeff(3,i)*t(j)^3 + (1/4)*coeff(4,i)*t(j)^4 - coeff(5,i)/t(j) + coeff(6,i) - coeff(7,i))... - -(coeff0(1,i)*t0 + (1/2)*coeff0(2,i)*t0^2 + (1/3)*coeff0(3,i)*t0^3 + (1/4)*coeff0(4,i)*t0^4 - coeff0(5,i)/t0 + coeff0(6,i) - coeff0(7,i)))... - * 1000; - end - end - - h_fc = Nf(4)*y_hf(4)+Nf(5)*y_hf(5)-N0(2)*y_hf(2); - del_hN = Nf(1).*del_h(1,:)+Nf(4).*del_h(4,:)+Nf(5).*del_h(5,:)+Nf(2).*del_h(2,:); - RT0 = sum(N0)*R*temp_0; - RTf = -sum(Nf)*R*Tf_guess; - - Tf_solve = h_fc + del_hN + RT0 + RTf'; - tol0 = Tf_guess(1) - Tf_guess(3); - - if Tf_solve(2) < 0 - if Tf_solve(1) > 0 - mp = 0.5*(Tf_guess(1)+Tf_guess(2)); - Tf_guess = [Tf_guess(1);mp;Tf_guess(2)]; - else - mp = 0.5*(Tf_guess(2)+Tf_guess(3)); - Tf_guess = [Tf_guess(2);mp;Tf_guess(3)]; - end - else - if Tf_solve(1) < 0 - mp = 0.5*(Tf_guess(1)+Tf_guess(2)); - Tf_guess = [Tf_guess(1);mp;Tf_guess(2)]; - else - mp = 0.5*(Tf_guess(2)+Tf_guess(3)); - Tf_guess = [Tf_guess(2);mp;Tf_guess(3)]; - end - end - - count = count+1; -end - -Tf=Tf_guess(2); -Pf=(Tf_guess(2)*sum(Nf)*R)/volume; -dP = Pf-pres_0; -Tf = Tf - 273.15; - -%----------------------- -% Create SS Vector of T and P -%----------------------- -tss=5*(1:16)+25; -Tf=Tf*ones(1,16); -dP=dP*ones(1,16); - -tss=[35;40;45;50;55;60]; -Tf=[Tf;Tf;Tf;Tf;Tf;Tf]; -dP=[dP;dP;dP;dP;dP;dP]; - -%------------------------------------------------ -%Write Expected Data CSV File -%------------------------------------------------ -yff(:,1) = tspan; - -header1 = {'Time','CH4','O2','CO2','H2O'}; -filename1 = '../../Verification/Species/reactionrate_EDC_1step_CH4_nonmix_spec.csv'; -fid = fopen(filename1,'wt'); -fprintf(fid,'%s, %s, %s, %s, %s\n',header1{:}); -for j=8:length(tspan) - fprintf(fid,'%f, %f, %f, %f, %f\n',yff(j,:)); -end -fclose(fid); - -header1 = {'Time','TEMP','PRES'}; -filename1 = '../../Verification/Species/reactionrate_EDC_1step_CH4_nonmix_temppres.csv'; -fid = fopen(filename1,'wt'); -fprintf(fid,'%s, %s, %s\n',header1{:}); -for j=1:length(tss) - fprintf(fid,'%f, %f, %f\n',tss(j),Tf(j),dP(j)); -end -fclose(fid); - - diff --git a/Utilities/Matlab/scripts/EDC_propane.m b/Utilities/Matlab/scripts/EDC_propane.m deleted file mode 100644 index 0d7d55d04be..00000000000 --- a/Utilities/Matlab/scripts/EDC_propane.m +++ /dev/null @@ -1,203 +0,0 @@ -%----------------------- -% C Weinschenk -% Verification of EDM (Species,Temp,Pres) -% Fuel: Propane -% 1/2012 -%----------------------- - -clear all -close all - -%----------------------- -% Initialize species -%----------------------- - -%------------- -% vector key -% (1) = nitrogen -% (2) = propane -% (3) = oxygen -% (4) = carbon dioxide -% (5) = water vapor -%------------- - -temp_0 = 298.15; % intial temperature [K] -volume = 0.001; % volume [m^3] -rho=1.2043424; % density [kg/m^3] -R = 8.3145; % gas constnat [J/mol K] -mass = rho* volume; % mass [kg] - -yf0 = [0.720828; 0.060321; 0.218851; 0.0; 0.0]; % initial mass fraction - -y_MW = [28.0134; 44.09562; 31.9988; 44.0095; 18.01528]; % [g/mol] - -y_hf = [0.0; -103850; 0.0; -393522; -241826]; % [J/mol] - -N0 = (1000*volume*rho*yf0)./(y_MW); % initial moles - -pres_0 = temp_0*sum(N0)*R/volume; % initial pressure - -%----------------------- -% Setup time vector for integration -%----------------------- -tspan=0:0.1:2; - -%----------------------- -% Setup vector of moles for integrator -%----------------------- -y0=[N0(1) N0(2) N0(3) N0(4) N0(5)]; - -%----------------------- -% Pass information to integrator -%----------------------- -options = odeset('RelTol',1e-10,'AbsTol',1e-10); -[T,Y]=ode113(@EDC_spec_propane,tspan,y0, options); - -Nf = Y(end,:); - -%----------------------- -% Convert back to mass fraction -%----------------------- -yff(:,1)=(y_MW(1)*Y(:,1))/(1000*rho*volume); -yff(:,2)=(y_MW(2)*Y(:,2))/(1000*rho*volume); -yff(:,3)=(y_MW(3)*Y(:,3))/(1000*rho*volume); -yff(:,4)=(y_MW(4)*Y(:,4))/(1000*rho*volume); -yff(:,5)=(y_MW(5)*Y(:,5))/(1000*rho*volume); - -%----------------------- -% Determine final temperature and pressure -%----------------------- -tol0 = 1e5; -tol = 0.1; -Tf_guess = [1000;2000;3000]; % final temperature guess [K] -count = 0; -while abs(tol0)>tol - %--------- - %NOTE: coefficients were found from NIST Webbook - %--------- - - %initial temperature coeffs - coeff0(:,1) = [28.98641;1.853978;-9.647459;16.63537;0.000117;-8.671914;0.0]; - coeff0(:,2) = [1;1;1;1;1;1;1]; %placeholder for propane - coeff0(:,3) = [31.32234;-20.23531;57.86644;-36.50624;-0.007374;-8.903471;0.0]; - coeff0(:,4) = [24.99735;55.18696;-33.69137;7.948387;-0.136638;-403.6075;-393.5224]; - coeff0(:,5) = [30.09200;6.832514;6.793435;-2.534480;0.082139;-250.8810;-241.8264]; - - %nitrogen cp coeffs [J/mol K] - if Tf_guess <=500 - coeff(:,1) = [28.98641;1.853978;-9.647459;16.63537;0.000117;-8.671914;0.0]; - elseif 500 < Tf_guess <= 2000 - coeff(:,1) = [19.50583;19.88705;-8.598535;1.369784;0.527601;-4.935202;0.0]; - else - coeff(:,1) = [35.51872;1.128728;-0.196103;0.014662;-4.553760;-18.97091;0.0]; - end - %propane cp coeffs [J/mol K] - if Tf_guess <= 1300 - coeff(:,2) = [1;1;1;1;1;1;1]; %placeholder for propane - else - coeff(:,2) = [1;1;1;1;1;1;1]; %placeholder for propane - end - %oxygen cp coeffs [J/mol K] - if Tf_guess <=700 - coeff(:,3) = [31.32234;-20.23531;57.86644;-36.50624;-0.007374;-8.903471;0.0]; - elseif 700 < Tf_guess <= 2000 - coeff(:,3) = [30.03235;8.772972;-3.988133;0.788313;-0.741599;-11.32468;0.0]; - else - coeff(:,3) = [20.91111;10.72071;-2.020498;0.146449;9.245722;5.337651;0.0]; - end - %carbon dooxide cp coeffs [J/mol K] - if Tf_guess <= 1200 - coeff(:,4) = [24.99735;55.18696;-33.69137;7.948387;-0.136638;-403.6075;-393.5224]; - else - coeff(:,4) = [58.16639;2.720074;-0.492289;0.038844;-6.447293;-425.9186;-393.5224]; - end - %water vapor cp coeffs [J/mol K] - if Tf_guess <= 1700 - coeff(:,5) = [30.09200;6.832514;6.793435;-2.534480;0.082139;-250.8810;-241.8264]; - else - coeff(:,5) = [41.96246;8.622053;-1.499780;0.098199;-11.15764;-272.1797;-241.8264]; - end - - t=Tf_guess./1000; - t0=temp_0/1000; - for j=1:3 - for i=1:5 - del_h(i,j) = ((coeff(1,i)*t(j) + (1/2)*coeff(2,i)*t(j)^2 + (1/3)*coeff(3,i)*t(j)^3 + (1/4)*coeff(4,i)*t(j)^4 - coeff(5,i)/t(j) + coeff(6,i) - coeff(7,i))... - -(coeff0(1,i)*t0 + (1/2)*coeff0(2,i)*t0^2 + (1/3)*coeff0(3,i)*t0^3 + (1/4)*coeff0(4,i)*t0^4 - coeff0(5,i)/t0 + coeff0(6,i) - coeff0(7,i)))... - * 1000; - end - end - - %overwrite del_h(2,:) for propane - coeff2(:,1) = [1.67536e-9;-5.46675e-6;4.38029e-3;2.7949;6.11728e2]; - for j=1:3 - del_h(2,j)=((1/5)*coeff2(1,1)*Tf_guess(j)^5+(1/4)*coeff2(2,1)*Tf_guess(j)^4+(1/3)*coeff2(3,1)*Tf_guess(j)^3+(1/2)*coeff2(4,1)*Tf_guess(j)^2+coeff2(5,1)*Tf_guess(j))... - -((1/5)*coeff2(1,1)*temp_0^5+(1/4)*coeff2(2,1)*temp_0^4+(1/3)*coeff2(3,1)*temp_0^3+(1/2)*coeff2(4,1)*temp_0^2+coeff2(5,1)*temp_0)... - .* (Tf_guess(j) - temp_0); - end - - h_fc = Nf(4)*y_hf(4)+Nf(5)*y_hf(5)-N0(2)*y_hf(2); - del_hN = Nf(1).*del_h(1,:)+Nf(4).*del_h(4,:)+Nf(5).*del_h(5,:); - RT0 = sum(N0)*R*temp_0; - RTf = -sum(Nf)*R*Tf_guess; - - Tf_solve = h_fc + del_hN + RT0 + RTf'; - tol0 = Tf_guess(1) - Tf_guess(3); - - if Tf_solve(2) < 0 - if Tf_solve(1) > 0 - mp = 0.5*(Tf_guess(1)+Tf_guess(2)); - Tf_guess = [Tf_guess(1);mp;Tf_guess(2)]; - else - mp = 0.5*(Tf_guess(2)+Tf_guess(3)); - Tf_guess = [Tf_guess(2);mp;Tf_guess(3)]; - end - else - if Tf_solve(1) < 0 - mp = 0.5*(Tf_guess(1)+Tf_guess(2)); - Tf_guess = [Tf_guess(1);mp;Tf_guess(2)]; - else - mp = 0.5*(Tf_guess(2)+Tf_guess(3)); - Tf_guess = [Tf_guess(2);mp;Tf_guess(3)]; - end - end - - count = count+1; -end -Tf=Tf_guess(2); -Pf=(Tf_guess(2)*sum(Nf)*R)/volume; -dP = Pf-pres_0; -Tf = Tf - 273.15; - -%----------------------- -% Create SS Vector of T and P -%----------------------- - -tss=[1;1.25;1.5;1.75;2]; -Tf=[Tf;Tf;Tf;Tf;Tf]; -dP=[dP;dP;dP;dP;dP]; - -%------------------------------------------------ -% Write Expected Data CSV File -%------------------------------------------------ -yff(:,1) = tspan; - -header1 = {'Time','O2','C3H8','CO2','H2O'}; -filename1 = '../../Verification/Species/reactionrate_EDC_flim_1step_C3H8_spec.csv'; -fid = fopen(filename1,'wt'); -fprintf(fid,'%s, %s, %s, %s, %s\n',header1{:}); -for j=1:length(tspan) - fprintf(fid,'%f, %f, %f, %f, %f\n',yff(j,:)); -end -fclose(fid); - -header1 = {'Time','TEMP','PRES'}; -filename1 = '../../Verification/Species/reactionrate_EDC_flim_1step_C3H8_temppres.csv'; -fid = fopen(filename1,'wt'); -fprintf(fid,'%s, %s, %s\n',header1{:}); -for j=1:length(tss) - fprintf(fid,'%f, %f, %f\n',tss(j),Tf(j),dP(j)); -end -fclose(fid); - - diff --git a/Utilities/Matlab/scripts/EDC_spec_methane.m b/Utilities/Matlab/scripts/EDC_spec_methane.m deleted file mode 100644 index 2a71a04f7fd..00000000000 --- a/Utilities/Matlab/scripts/EDC_spec_methane.m +++ /dev/null @@ -1,35 +0,0 @@ -function dy=EDC_spec_methane(t,y) - -% (1) - nitrogen -% (2) - methane -% (3) - oxygen -% (4) - co2 -% (5) - h2o - -%----------------------- -% Initialize variables -%----------------------- - -m=1; % chemical reaction coefficient fuel -s=2; % chemical reaction coefficient oxygen -c=1; % chemical reaction coefficient carbon dioxide -h=2; % chemical reaction coefficient water vapor -tau_mix=0.125; % set mixing time [s] - -%----------------------- -% Logic for setting molar production rate -% No suppression in this case -%----------------------- -r1=min([y(2) y(3)/s]); -r_co2(1)=m*r1/tau_mix; - -%----------------------- -% Integration -%----------------------- -dy(1)= 0.0; -dy(2)=-(1/m)*r_co2(1); -dy(3)=-(s/m)*r_co2(1); -dy(4)= c*r_co2(1); -dy(5)= h*r_co2(1); - -dy=[dy(1);dy(2);dy(3);dy(4);dy(5)]; \ No newline at end of file diff --git a/Utilities/Matlab/scripts/EDC_spec_propane.m b/Utilities/Matlab/scripts/EDC_spec_propane.m deleted file mode 100644 index d0c65643cb9..00000000000 --- a/Utilities/Matlab/scripts/EDC_spec_propane.m +++ /dev/null @@ -1,35 +0,0 @@ -function dy=EDC_spec_propane(t,y) - -% (1) - nitrogen -% (2) - propane -% (3) - oxygen -% (4) - co2 -% (5) - h2o - -%----------------------- -% Initialize variables -%----------------------- - -m=1; % chemical reaction coefficient fuel -s=5; % chemical reaction coefficient oxygen -c=3; % chemical reaction coefficient carbon dioxide -h=4; % chemical reaction coefficient water vapor -tau_mix=0.125; % set mixing time [s] - -%----------------------- -% Logic for setting molar production rate -% No suppression in this case -%----------------------- -r1=min([y(2) y(3)/s]); -r_co2(1)=m*r1/tau_mix; - -%----------------------- -% Integration -%----------------------- -dy(1)= 0.0; -dy(2)=-(1/m)*r_co2(1); -dy(3)=-(s/m)*r_co2(1); -dy(4)= c*r_co2(1); -dy(5)= h*r_co2(1); - -dy=[dy(1);dy(2);dy(3);dy(4);dy(5)]; \ No newline at end of file diff --git a/Utilities/Matlab/scripts/EDC_species.m b/Utilities/Matlab/scripts/EDC_species.m deleted file mode 100644 index 183468068bc..00000000000 --- a/Utilities/Matlab/scripts/EDC_species.m +++ /dev/null @@ -1,12 +0,0 @@ -%----------------------- -% C Weinschenk -% Calls scripts for EDM verfication -% 1/2012 -%----------------------- - -close all -clear all - -EDC_methane; % calls EDM_spec_methane.m function file for species integration -EDC_propane; % calls EDM_spec_propane.m function file for species integration -EDC_methane_nonmix; % calls EDM_spec_methane.m function file for species integration \ No newline at end of file diff --git a/Utilities/Matlab/scripts/Monin_Obukhov_Profiles.m b/Utilities/Matlab/scripts/Monin_Obukhov_Profiles.m deleted file mode 100644 index 6a0e7e717e1..00000000000 --- a/Utilities/Matlab/scripts/Monin_Obukhov_Profiles.m +++ /dev/null @@ -1,104 +0,0 @@ -% McGrattan -% 08-17-2017 -% Monin_Obukhov_Profiles.m -% -% Atmospheric Boundary Layer profiles, based on M-O theory as described in -% "Falcon Series Data Report", GRI-89/0138, June 1990. -% These plots are used as illustrations in the FDS User's Guide. - -clear all -close all - -plot_style - -rho_0 = 1.2; -g = 9.81; -cp = 1005; -kappa = 0.41; -p_0 = 100000; - -L = -100.; -z_0 = 0.01; -z_r = 2.; -u_r = 5.; -p_r = p_0-rho_0*g*(z_r-z_0); -T_r = 20+273.15; -theta_r = T_r*(p_0/p_r)^0.285; -if L<0; - x_r = (1-16*z_r/L).^0.25; - psi_m_r = 2*log((1+x_r)/2) + log((1+x_r^2)/2) - 2*atan(x_r) + pi/2; - psi_h_r = 2*log((1+x_r^2)/2); -else - psi_m_r = -5*z_r/L; - psi_h_r = psi_m_r; -end -u_star = kappa*u_r/(log(z_r/z_0)-psi_m_r); -theta_0 = theta_r/(1+u_star^2*(log(z_r/z_0)-psi_h_r)/(g*kappa^2*L)); -theta_star = u_star^2*theta_0/(g*kappa*L); - -z = [z_0 10*z_0 1 2 3 4 5 6 7 8 9 10 15 20 25 30 50 100]; - -figure(1) -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - -if L<0; -x = (1-16*z/L).^0.25; -psi_m = 2*log((1+x)/2) + log((1+x.^2)/2) - 2*atan(x) + pi/2; -psi_h = 2*log((1+x.^2)/2); -else -psi_m = -5*z/L; -psi_h = psi_m; -end - -u = (u_star/kappa)*(log(z/z_0) - psi_m); -theta = theta_0 + (theta_star/kappa)*(log(z/z_0) - psi_h); -T = theta.*(p_0./(p_0-rho_0*g*(z-z_0))).^-0.285; - -plot(u,z,'ko-'); hold on -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Title_Font_Size) -axis([0 10 0 30]) -text(.05,.90,'$u(2 \; \hbox{m})=5$ m/s','FontName',Font_Name,'FontSize',Label_Font_Size,'Interpreter','LaTeX','Units','normalized') -text(.05,.80,['$L=' num2str(L,'%4.0f\n') '$ m'],'FontName',Font_Name,'FontSize',Label_Font_Size,'Interpreter','LaTeX','Units','normalized') -text(.05,.70,['$z_0=' num2str(z_0,'%5.3f\n') '$ m'],'FontName',Font_Name,'FontSize',Label_Font_Size,'Interpreter','LaTeX','Units','normalized') -xlabel('Velocity (m/s)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter) -ylabel('Height (m)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter) - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); -print(gcf,'-dpdf',['vel_' 'L=' num2str(L,'%4.0f\n')]) - -figure(2) -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - -plot(T-273.15,z,'ko-'); hold on -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Title_Font_Size) -axis([15 25 0 30]) -text(.05,.90,'$T(2 \; \hbox{m})=20\;^\circ$C','FontName',Font_Name,'FontSize',Label_Font_Size,'Interpreter','LaTeX','Units','normalized') -text(.05,.80,['$L=' num2str(L,'%4.0f\n') '$ m'],'FontName',Font_Name,'FontSize',Label_Font_Size,'Interpreter','LaTeX','Units','normalized') -text(.05,.70,['$z_0=' num2str(z_0,'%5.3f\n') '$ m'],'FontName',Font_Name,'FontSize',Label_Font_Size,'Interpreter','LaTeX','Units','normalized') -xlabel('Temperature (°C)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter) -ylabel('Height (m)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter) - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); -print(gcf,'-dpdf',['tmp_' 'L=' num2str(L,'%4.0f\n')]) - -% Print out the u and T profiles - -fid = fopen('Monin_Obukhov_Profile.csv','wt','n'); -fprintf(fid,'%s\n','z,u,T'); -for j=1:numel(z) - fprintf(fid,'%6.2f, %5.2f, %6.2f\n',z(j),u(j),T(j)-273.15); -end -fclose(fid); - diff --git a/Utilities/Matlab/scripts/backward_facing_step.m b/Utilities/Matlab/scripts/backward_facing_step.m deleted file mode 100644 index e6cfcc28f28..00000000000 --- a/Utilities/Matlab/scripts/backward_facing_step.m +++ /dev/null @@ -1,892 +0,0 @@ -% Toms -% 8-8-14 -% backward_facing_step.m - -function main() - -clear all -close all - -disp('backward_facing_step ...') - -expdir = '../../../exp/Backward_Facing_Step/'; -datdir = '../../../out/Backward_Facing_Step/'; -pltdir = '../../Manuals/FDS_Validation_Guide/SCRIPT_FIGURES/Backward_Facing_Step/'; - -rkappa = 1/.41; -B = 5.2; -rho = 1.1992662; % density, kg/m^3 -U_0 = 7.72; % reference velocity, m/s -h = 0.0098; % step height, m -dx = h/20; -visc = 1.7801E-5; -visc_nu = visc/rho; -x_loc = {'-3','4','6','10'}; - -nx = [5 10 20]; -lnx = length(nx); -dx = h./nx; -fds_marker = {'b-' 'r-' 'm-'}; -fds_key = {'{\it h/\deltax}=5' '{\it h/\deltax}=10' '{\it h/\deltax}=20'}; - -if ~exist([expdir,'backward_facing_step_data.csv']) - display(['Error: File ' [expdir,'backward_facing_step_data.csv'] ' does not exist. Skipping case.']) - return -end -for i=1:lnx - if ~exist([datdir,'backward_facing_step_',num2str(nx(i)),'_line.csv']) - display(['Error: File ' [datdir,'backward_facing_step_',num2str(nx(i)),'_line.csv'] ' does not exist. Skipping case.']) - return - end -end - -% read experimental and FDS data files - -D = importdata([expdir,'backward_facing_step_data.csv'],',',1); -for i=1:lnx - M{i} = importdata([datdir,'backward_facing_step_',num2str(nx(i)),'_line.csv'],',',2); -end - -linetype_vect_leg = [{'-'},{'-'},{'-'}]; -linetype_vect_noleg = [{'-'},{'-'},{'-'}]; -symbol_vect_noleg = [{'p'},{'s'},{'d'}]; - -line_color_vect={[20/255,168/255,113/225],[186/255,62/255,62/255],[64/255,97/255,191/255]}; - -plot_style - -f1 = figure(1); -set(f1,'Visible',Figure_Visibility) -p1 = tight_subplot(1,1, [.01 .01],[.13 .08],[.11 .019]); - -f2 = figure(2); -set(f2,'Visible',Figure_Visibility); -p2 = tight_subplot(1,1, [.01 .01],[.13 .08],[.141 .019]); - -f3 = figure(3); -set(f3,'Visible',Figure_Visibility); -sp1 = tight_subplot(1,4, [.01 .01],[.142 .055],[.108 .01]); - -f4 = figure(4); -set(f4,'Visible',Figure_Visibility); -sp2 = tight_subplot(1,4, [.01 .01],[.142 .055],[.108 .01]); - -f5 = figure(5); -set(f5,'Visible',Figure_Visibility); -sp3 = tight_subplot(1,4, [.01 .01],[.142 .055],[.108 .01]); - -f6 = figure(6); -set(f6,'Visible',Figure_Visibility); -sp4 = tight_subplot(1,4, [.01 .01],[.142 .055],[.108 .01]); - -f7 = figure(7); -set(f7,'Visible',Figure_Visibility); -sp5 = tight_subplot(1,4, [.01 .01],[.142 .055],[.108 .01]); - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% streamwise data along bottom of channel for determining reattachment location -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -set(0, 'CurrentFigure', f1) - -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - -%%%Exp Data%%% - -j = find(strcmp(D.colheaders,'Cf-x/h')); -xoh = D.data(:,j); -j = find(strcmp(D.colheaders,'Cf')); -Cf = D.data(:,j); - -error = .0005 * ones(length(Cf),1); - -h_dat=errorbar(xoh,Cf,error,'ko','markersize',10); - -hold on - -axis([0 20 -4E-3 5E-3]) -set(gca,'YTick',[-4E-3:1E-3:5E-3]) -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Label_Font_Size) - -%%%FDS Data%%% - -for q = 1:lnx - j= find(strcmp(M{q}.colheaders, 'u_tau-x')); - x = M{q}.data(:,j); - I = find(x > 0); - x = x(I); - - j = find(strcmp(M{q}.colheaders, 'u_tau')); - u_tau = M{q}.data(I,j); - tau_w = rho*u_tau.^2; - - j = find(strcmp(M{q}.colheaders, 'u_wall')); - u_wall = M{q}.data(I,j); - - Cf_fds = visc*u_wall/(.5*h/nx(q))/(.5*rho*U_0^2); - - plot(x/h,Cf_fds, char(linetype_vect_leg(q)), 'Color',line_color_vect{q}, 'LineWidth', 1); - - h_leg(q) = plot(x(1:round(length(x)/15):end)/h, Cf_fds(1:round(length(Cf_fds)/15):end), char(symbol_vect_noleg(q)), 'Color', line_color_vect{q}, 'LineWidth', 1, 'MarkerSize', 10); - - yh = ylabel('{\it C_f}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - xh = xlabel('{\it x/h}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - - % if q == 1 - % set(xh,'Units','Pixels') - % xh_pos = get(xh,'Position'); - % xh_pos(1) = xh_pos(1)*1.073; - % xh_pos(2) = xh_pos(2)*1.17; - % set(xh,'Position',xh_pos) - % end - - if q == 1 - Pos=get(yh,'Position'); - set(yh,'Position',[Pos(1)*.95 Pos(2)*.5 Pos(3)]) - end -end - -axis([0 20 -4E-3 4E-3]) - -lh = legend([h_dat,h_leg([1:length(h_leg)])], ['J&D',fds_key([1:length(h_leg)])], 'Location', 'SouthEast'); -set(lh,'box','off') -%set(lh,'Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - -% add Git revision if file is available - -Git_Filename = [datdir,'backward_facing_step_',num2str(nx(1)),'_git.txt']; -addverstr(gca,Git_Filename,'linear') - -% print to pdf - -print(gcf,'-dpdf',[pltdir,'backward_facing_step_Cf']) - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% pressure coefficient along bottom of channel -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -set(0, 'CurrentFigure', f2) - -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X,Plot_Y,Plot_Width,Plot_Height]) -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - -%%%Exp Data%%% -j = find(strcmp(D.colheaders,'Cp-x/h')); -xoh = D.data(:,j); -I = find(xoh<=20); -xoh = xoh(I); -j = find(strcmp(D.colheaders,'Cp')); -Cp = D.data(I,j); - -h_dat=plot(xoh,Cp,'ko','markersize',10); - -hold on - -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Label_Font_Size) - -%%%FDS Data%%% - -for q = 1:lnx - j = find(strcmp(M{q}.colheaders, 'cp-x')); - x = M{q}.data(:,j); - I = find(x>0); - x = x(I); - - j = find(strcmp(M{q}.colheaders, 'cp')); - Cp_fds = M{q}.data(I,j)/(U_0^2); - Cp_fds = Cp_fds + Cp(end) - Cp_fds(end); - - plot(x/h,Cp_fds, char(linetype_vect_leg(q)), 'Color',line_color_vect{q}, 'LineWidth', 1); - - h_leg2(q) = plot(x(1:round(length(x)/15):end)/h, Cp_fds(1:round(length(Cp_fds)/15):end), char(symbol_vect_noleg(q)), 'Color', line_color_vect{q}, 'LineWidth', 1, 'MarkerSize', 10); - - yh = ylabel('{\it C_p}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - xh = xlabel('\it x/h','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - - Pos=get(yh, 'Position'); - set(yh, 'Position',[Pos(1)*.93 Pos(2) Pos(3)]); -end - -lh = legend([h_dat,h_leg2([1:length(h_leg2)])], ['J&D',fds_key([1:length(h_leg2)])], 'Location', 'SouthEast'); - -% add version string if file is available - -Git_Filename = [datdir,'backward_facing_step_',num2str(nx(1)),'_git.txt']; -addverstr(gca,Git_Filename,'linear') - -% print to pdf - -print(gcf,'-dpdf',[pltdir,'backward_facing_step_Cp']) - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% wall-normal streamwise velocity profiles -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -set(0, 'CurrentFigure', f3); - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - -for i = 1:4 - - %%%Exp Data%%% - - j = find(strcmp(D.colheaders,strcat({'z '},x_loc{i}))); - z = D.data(:,j); - I = find(z>0); - if i == 1 - z=z(I)*0.001+h; - else - z=z(I)*0.001; - end - j = find(strcmp(D.colheaders,strcat({'U '},x_loc{i}))); - u_data = D.data(I,j); - - set(f3,'CurrentAxes',sp1(i)) - H(1)=plot(u_data/U_0,z/h,'ko','markersize',10); - hold on - - %%%FDS data%%% - - for q=1:lnx - j = find(strcmp(M{q}.colheaders, strcat({'U-VEL '},x_loc{i},'-z'))); - z = M{q}.data(:,j); - I = find(z>0); %ensures no data for z<0 is incorporated... - z = z(I); - - j = find(strcmp(M{q}.colheaders, strcat({'U-VEL '},x_loc{i}))); - u_vel = M{q}.data(I,j); - - plot(u_vel/U_0,z/h,char(linetype_vect_leg(q)), 'Color',line_color_vect{q}, 'LineWidth', 1) - plot(u_vel(1:round(length(u_vel)/15):end)/U_0,z(1:round(length(z)/15):end)/h,char(symbol_vect_noleg(q)), 'MarkerSize', 10, 'Color',line_color_vect{q}, 'LineWidth', 1) - - axis([-.2 1.0 0 3.5]) - - set(gca,'FontName',Font_Name) - set(gca,'FontSize',Label_Font_Size) - - Pos = get(gca,'Position'); - set(gca,'Position',[Pos(1) Pos(2)+.015 Pos(3) Pos(4)*0.98]) - - if i == 1 - ylabel('{\it z/h}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - xh=xlabel('{\it /U_0}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - else - set(gca,'XTickLabel',[]) - set(gca,'YTickLabel',[]) - end - if q == 1 - tpos = [0.01 3.2]; % from get(th,'Position') where th is a title handle - th = text(tpos(1),tpos(2),strcat('{\itx/h}','=',x_loc{i}),'Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - set(th,'Position',tpos) - end - end - -end - -% add Git revision if file is available - -Git_Filename = [datdir,'backward_facing_step_',num2str(nx(1)),'_git.txt']; -addverstr(gca,Git_Filename,'linear',0,1.05) - -% print to pdf - -print(gcf,'-dpdf',[pltdir,'backward_facing_step_U']) - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% wall-normal vertical velocity profiles -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -set(0, 'CurrentFigure', f4); - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - -for i = 1:4 - - j = find(strcmp(D.colheaders,strcat({'z '},x_loc{i}))); - z = D.data(:,j); - I = find(z>0); - if i == 1 - z=z(I)*0.001+h; - else - z=z(I)*0.001; - end - j = find(strcmp(D.colheaders,strcat({'W '},x_loc{i}))); - w_data = D.data(I,j); - - set(f4,'CurrentAxes',sp2(i)) - plot(w_data/U_0,z/h,'ko','markersize',10); - hold on - - for q=1:lnx - j = find(strcmp(M{q}.colheaders,strcat({'W-VEL '},x_loc{i},'-z'))); - z = M{q}.data(:,j); - I = find(z>0); - z = z(I); - - j = find(strcmp(M{q}.colheaders,strcat({'W-VEL '},x_loc{i}))); - w_vel = M{q}.data(I,j); - - plot(w_vel/U_0,z/h,char(linetype_vect_noleg(q)), 'MarkerSize', 1.5,'Color',line_color_vect{q}, 'LineWidth', 1); - plot(w_vel(1:round(length(w_vel)/15):end)/U_0,z(1:round(length(z)/15):end)/h, char(symbol_vect_noleg(q)), 'MarkerSize', 10, 'Color',line_color_vect{q}, 'LineWidth', 1) - - axis([-.1 .1 0 3.5]) - - set(gca,'FontName',Font_Name) - set(gca,'FontSize',Label_Font_Size) - - Pos = get(gca,'Position'); - set(gca,'Position',[Pos(1) Pos(2)+.015 Pos(3) Pos(4)*0.98]) - - if i == 1 - ylabel('{\it z/h}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - xh=xlabel('{\it /U_0}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - else - set(gca,'XTickLabel',[]) - set(gca,'YTickLabel',[]) - end - - if q == 1 - tpos = [0.01 3.2]; % from get(th,'Position') where th is a title handle - th = text(tpos(1),tpos(2),strcat('{\itx/h}','=',x_loc{i}),'Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - set(th,'Position',tpos) - end - end - - %th = title(strcat('{\it x/h}','=',x_loc{i}),'Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name, 'HorizontalAlignment', 'Left'); - -end - -% add Git revision if file is available - -Git_Filename = [datdir,'backward_facing_step_',num2str(nx(1)),'_git.txt']; -addverstr(gca,Git_Filename,'linear',0,1.05) - -% print to pdf - -print(gcf,'-dpdf',[pltdir,'backward_facing_step_W']) - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% wall-normal uu profiles -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -set(0, 'CurrentFigure', f5); - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - -for i = 1:4 - - j = find(strcmp(D.colheaders,strcat({'z '}, x_loc{i}))); - z = D.data(:,j); - I = find(z>0); - if i == 1 - z=z(I)*0.001+h; - else - z=z(I)*0.001; - end - j = find(strcmp(D.colheaders,strcat({'uu '},x_loc{i}))); - uu_data = D.data(I,j); - - error = uu_data * .15 / U_0^2; - - set(f5,'CurrentAxes',sp3(i)) - plot(uu_data/U_0^2,z/h,'ko','markersize',10); - hold on - %herrorbar(uu_data/U_0^2,z/h, error,'ko'); - - for q=1:lnx - - j = find(strcmp(M{q}.colheaders,strcat({'uu '},x_loc{i},'-z'))); - z = M{q}.data(:,j); - I = find(z>0); - z = z(I); - - j = find(strcmp(M{q}.colheaders,strcat({'uu '},x_loc{i}))); - uu_fds = M{q}.data(I,j); - - plot(uu_fds/U_0^2,z/h, char(linetype_vect_noleg(q)), 'MarkerSize', 2, 'Color',line_color_vect{q}, 'LineWidth', 1); - plot(uu_fds(1:round(length(uu_fds)/15):end)/U_0^2,z(1:round(length(z)/15):end)/h,char(symbol_vect_noleg(q)), 'MarkerSize', 10, 'Color',line_color_vect{q}, 'LineWidth', 1) - - axis([0 0.04 0 3.5]) - - set(gca,'FontName',Font_Name) - set(gca,'FontSize',Label_Font_Size) - - Pos = get(gca,'Position'); - set(gca,'Position',[Pos(1) Pos(2)+.015 Pos(3) Pos(4)*0.98]) - - if i == 1 - ylabel('{\it z/h}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - xh=xlabel('{\it /U_0^2}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - else - set(gca,'YTickLabel',[]) - set(gca,'XTickLabel',[]) - end - if q == 1 - tpos = [0.01 3.2]; % from get(th,'Position') where th is a title handle - th = text(tpos(1),tpos(2),strcat('{\itx/h}','=',x_loc{i}),'Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - set(th,'Position',tpos) - end - end - -end - -% add Git revision if file is available - -Git_Filename = [datdir,'backward_facing_step_',num2str(nx(1)),'_git.txt']; -addverstr(gca,Git_Filename,'linear',0,1.05) - -% print to pdf - -print(gcf,'-dpdf',[pltdir,'backward_facing_step_uu']) - - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% wall-normal ww profiles -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -set(0, 'CurrentFigure', f6); - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - -for i = 1:4 - - j = find(strcmp(D.colheaders,strcat({'z '},x_loc{i}))); - z = D.data(:,j); - I = find(z>0); - if i == 1 - z=z(I)*0.001+h; - else - z=z(I)*0.001; - end - j = find(strcmp(D.colheaders,strcat({'ww '},x_loc{i}))); - ww_data = D.data(I,j); - - error = ww_data * .15 / U_0^2; - - set(f6,'CurrentAxes',sp4(i)) - plot(ww_data/U_0^2,z/h,'ko','markersize',10); - hold on - %herrorbar(ww_data/U_0^2,z/h,error,'ko'); - - for q = 1:lnx - j = find(strcmp(M{q}.colheaders,strcat({'ww '},x_loc{i},'-z'))); - z = M{q}.data(:,j); - I = find(z>0); - z = z(I); - - - j = find(strcmp(M{q}.colheaders,strcat({'ww '},x_loc{i}))); - ww_fds = M{q}.data(I,j); - - - plot(ww_fds/U_0^2,z/h,char(linetype_vect_noleg(q)), 'MarkerSize', 1.5,'Color',line_color_vect{q}, 'LineWidth', 1); - plot(ww_fds(1:round(length(ww_fds)/15):end)/U_0^2,z(1:round(length(z)/15):end)/h,char(symbol_vect_noleg(q)), 'MarkerSize', 10, 'Color',line_color_vect{q}, 'LineWidth', 1) - - axis([0 .04 0 3.5]) - - set(gca,'FontName',Font_Name) - set(gca,'FontSize',Label_Font_Size) - - Pos = get(gca,'Position'); - set(gca,'Position',[Pos(1) Pos(2)+.015 Pos(3) Pos(4)*0.98]) - - if i == 1 - ylabel('{\it z/h}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name) - xh=xlabel('{\it /U_0^2}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - else - set(gca,'YTickLabel',[]) - set(gca,'XTickLabel',[]) - end - - if q == 1 - tpos = [0.01 3.2]; % from get(th,'Position') where th is a title handle - th = text(tpos(1),tpos(2),strcat('{\itx/h}','=',x_loc{i}),'Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - set(th,'Position',tpos) - end - end - -end - -% add Git if file is available - -Git_Filename = [datdir,'backward_facing_step_',num2str(nx(1)),'_git.txt']; -addverstr(gca,Git_Filename,'linear',0,1.05) - -% print to pdf - -print(gcf,'-dpdf',[pltdir,'backward_facing_step_ww']) - -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -% wall-normal uw profiles -%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - -set(0, 'CurrentFigure', f7); - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - -for i = 1:4 - j = find(strcmp(D.colheaders,strcat({'z '},x_loc{i}))); - z = D.data(:,j); - I = find(z>0); - if i == 1 - z=z(I)*0.001+h; - else - z=z(I)*0.001; - end - j = find(strcmp(D.colheaders,strcat({'uw '},x_loc{i}))); - uw_data = D.data(I,j); - - error = uw_data * .15 / U_0^2; - - set(f7,'CurrentAxes',sp5(i)) - h_dat=plot(uw_data/U_0^2,z/h,'ko','markersize',10); - hold on - %herrorbar(uw_data/U_0^2,z/h,error,'ko'); - - for q = 1:lnx - - j = find(strcmp(M{q}.colheaders,strcat({'uw '},x_loc{i},'-z'))); - z = M{q}.data(:,j); - I = find(z>0); - z = z(I); - - j = find(strcmp(M{q}.colheaders,strcat({'uw '},x_loc{i}))); - uw_fds = -1*M{q}.data(I,j); - - - plot(uw_fds/U_0^2,z/h,char(linetype_vect_noleg(q)), 'MarkerSize', 1.5,'Color',line_color_vect{q}, 'LineWidth', 1); - h_leg(q)=plot(uw_fds(1:round(length(uw_fds)/15):end)/U_0^2,z(1:round(length(z)/15):end)/h,char(symbol_vect_noleg(q)), 'MarkerSize', 10, 'Color',line_color_vect{q}, 'LineWidth', 1); - - axis([0 .04 0 3.5]) - - set(gca,'FontName',Font_Name) - set(gca,'FontSize',Label_Font_Size) - - Pos = get(gca,'Position'); - set(gca,'Position',[Pos(1) Pos(2)+.015 Pos(3) Pos(4)*0.98]) - - if i == 1 - ylabel('{\it z/h}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name) - xh=xlabel('{\it -/U_0^2}','Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - else - set(gca,'YTickLabel',[]) - set(gca,'XTickLabel',[]) - end - - if q == 1 - tpos = [0.01 3.2]; % from get(th,'Position') where th is a title handle - th = text(tpos(1),tpos(2),strcat('{\itx/h}','=',x_loc{i}),'Interpreter',Font_Interpreter,'Fontsize',Label_Font_Size,'FontName',Font_Name); - set(th,'Position',tpos) - end - end - -end - -% add version string if file is available - -Git_Filename = [datdir,'backward_facing_step_',num2str(nx(1)),'_git.txt']; -addverstr(gca,Git_Filename,'linear',0,1.05) - -% print to pdf - -print(gcf,'-dpdf',[pltdir,'backward_facing_step_uw']) - -end % main() - - -%THE FOLLOWING CODE IS ALL EXTERNAL SCRIPTS IMBEDDED IN THIS SCRIPT TO -%LIMIT THE TOMS INFLUENCE ON THE REPOSITORY - -function ha = tight_subplot(Nh, Nw, gap, marg_h, marg_w) - -% tight_subplot creates "subplot" axes with adjustable gaps and margins -% -% ha = tight_subplot(Nh, Nw, gap, marg_h, marg_w) -% -% in: Nh number of axes in hight (vertical direction) -% Nw number of axes in width (horizontaldirection) -% gap gaps between the axes in normalized units (0...1) -% or [gap_h gap_w] for different gaps in height and width -% marg_h margins in height in normalized units (0...1) -% or [lower upper] for different lower and upper margins -% marg_w margins in width in normalized units (0...1) -% or [left right] for different left and right margins -% -% out: ha array of handles of the axes objects -% starting from upper left corner, going row-wise as in -% going row-wise as in -% -% Example: ha = tight_subplot(3,2,[.01 .03],[.1 .01],[.01 .01]) -% for ii = 1:6; axes(ha(ii)); plot(randn(10,ii)); end -% set(ha(1:4),'XTickLabel',''); set(ha,'YTickLabel','') - -% Pekka Kumpulainen 20.6.2010 @tut.fi -% Tampere University of Technology / Automation Science and Engineering - - -if nargin<3; gap = .02; end -if nargin<4 || isempty(marg_h); marg_h = .05; end -if nargin<5; marg_w = .05; end - -if numel(gap)==1; - gap = [gap gap]; -end -if numel(marg_w)==1; - marg_w = [marg_w marg_w]; -end -if numel(marg_h)==1; - marg_h = [marg_h marg_h]; -end - -axh = (1-sum(marg_h)-(Nh-1)*gap(1))/Nh; -axw = (1-sum(marg_w)-(Nw-1)*gap(2))/Nw; - -py = 1-marg_h(2)-axh; - -ha = zeros(Nh*Nw,1); -ii = 0; -for ih = 1:Nh - px = marg_w(1); - - for ix = 1:Nw - ii = ii+1; - ha(ii) = axes('Units','normalized', ... - 'Position',[px py axw axh], ... - 'XTickLabel','', ... - 'YTickLabel',''); - px = px+axw+gap(2); - end - py = py-axh-gap(1); -end - -end - - -function hh = herrorbar(x, y, l, u, symbol) -%HERRORBAR Horizontal Error bar plot. -% HERRORBAR(X,Y,L,R) plots the graph of vector X vs. vector Y with -% horizontal error bars specified by the vectors L and R. L and R contain the -% left and right error ranges for each point in X. Each error bar -% is L(i) + R(i) long and is drawn a distance of L(i) to the right and R(i) -% to the right the points in (X,Y). The vectors X,Y,L and R must all be -% the same length. If X,Y,L and R are matrices then each column -% produces a separate line. -% -% HERRORBAR(X,Y,E) or HERRORBAR(Y,E) plots X with error bars [X-E X+E]. -% HERRORBAR(...,'LineSpec') uses the color and linestyle specified by -% the string 'LineSpec'. See PLOT for possibilities. -% -% H = HERRORBAR(...) returns a vector of line handles. -% -% Example: -% x = 1:10; -% y = sin(x); -% e = std(y)*ones(size(x)); -% herrorbar(x,y,e) -% draws symmetric horizontal error bars of unit standard deviation. -% -% This code is based on ERRORBAR provided in MATLAB. -% -% See also ERRORBAR - -% Jos van der Geest -% email: jos@jasen.nl -% -% File history: -% August 2006 (Jos): I have taken back ownership. I like to thank Greg Aloe from -% The MathWorks who originally introduced this piece of code to the -% Matlab File Exchange. -% September 2003 (Greg Aloe): This code was originally provided by Jos -% from the newsgroup comp.soft-sys.matlab: -% http://newsreader.mathworks.com/WebX?50@118.fdnxaJz9btF^1@.eea3ff9 -% After unsuccessfully attempting to contact the orignal author, I -% decided to take ownership so that others could benefit from finding it -% on the MATLAB Central File Exchange. - -if min(size(x))==1, - npt = length(x); - x = x(:); - y = y(:); - if nargin > 2, - if ~isstr(l), - l = l(:); - end - if nargin > 3 - if ~isstr(u) - u = u(:); - end - end - end -else - [npt,n] = size(x); -end - -if nargin == 3 - if ~isstr(l) - u = l; - symbol = '-'; - else - symbol = l; - l = y; - u = y; - y = x; - [m,n] = size(y); - x(:) = (1:npt)'*ones(1,n);; - end -end - -if nargin == 4 - if isstr(u), - symbol = u; - u = l; - else - symbol = '-'; - end -end - -if nargin == 2 - l = y; - u = y; - y = x; - [m,n] = size(y); - x(:) = (1:npt)'*ones(1,n);; - symbol = '-'; -end - -u = abs(u); -l = abs(l); - -if isstr(x) | isstr(y) | isstr(u) | isstr(l) - error('Arguments must be numeric.') -end - -if ~isequal(size(x),size(y)) | ~isequal(size(x),size(l)) | ~isequal(size(x),size(u)), - error('The sizes of X, Y, L and U must be the same.'); -end - -tee = (max(y(:))-min(y(:)))/100; % make tee .02 x-distance for error bars -% changed from errorbar.m -xl = x - l; -xr = x + u; -ytop = y + tee; -ybot = y - tee; -n = size(y,2); -% end change - -% Plot graph and bars -hold_state = ishold; -cax = newplot; -next = lower(get(cax,'NextPlot')); - -% build up nan-separated vector for bars -% changed from errorbar.m -xb = zeros(npt*9,n); -xb(1:9:end,:) = xl; -xb(2:9:end,:) = xl; -xb(3:9:end,:) = NaN; -xb(4:9:end,:) = xl; -xb(5:9:end,:) = xr; -xb(6:9:end,:) = NaN; -xb(7:9:end,:) = xr; -xb(8:9:end,:) = xr; -xb(9:9:end,:) = NaN; - -yb = zeros(npt*9,n); -yb(1:9:end,:) = ytop; -yb(2:9:end,:) = ybot; -yb(3:9:end,:) = NaN; -yb(4:9:end,:) = y; -yb(5:9:end,:) = y; -yb(6:9:end,:) = NaN; -yb(7:9:end,:) = ytop; -yb(8:9:end,:) = ybot; -yb(9:9:end,:) = NaN; -% end change - - -[ls,col,mark,msg] = colstyle(symbol); if ~isempty(msg), error(msg); end -symbol = [ls mark col]; % Use marker only on data part -esymbol = ['-' col]; % Make sure bars are solid - -h = plot(xb,yb,esymbol); hold on -h = [h;plot(x,y,symbol)]; - -if ~hold_state, hold off; end - -if nargout>0, hh = h; end - -end - -function cell2csv(filename,cellArray,delimiter) -% Writes cell array content into a *.csv file. -% -% CELL2CSV(filename,cellArray,delimiter) -% -% filename = Name of the file to save. [ i.e. 'text.csv' ] -% cellarray = Name of the Cell Array where the data is in -% delimiter = seperating sign, normally:',' (default) -% -% by Sylvain Fiedler, KA, 2004 -% modified by Rob Kohr, Rutgers, 2005 - changed to english and fixed delimiter -if nargin<3 - delimiter = ','; -end - -datei = fopen(filename,'w'); -for z=1:size(cellArray,1) - for s=1:size(cellArray,2) - - var = eval(['cellArray{z,s}']); - - if size(var,1) == 0 - var = ''; - end - - if isnumeric(var) == 1 - var = num2str(var); - end - - fprintf(datei,var); - - if s ~= size(cellArray,2) - fprintf(datei,[delimiter]); - end - end - fprintf(datei,'\n'); -end -fclose(datei); -end diff --git a/Utilities/Matlab/scripts/check_hrr.m b/Utilities/Matlab/scripts/check_hrr.m deleted file mode 100644 index 2498b7d6f8b..00000000000 --- a/Utilities/Matlab/scripts/check_hrr.m +++ /dev/null @@ -1,65 +0,0 @@ -% McDermott -% 9-20-10 -% check_hrr.m - -close all -clear all - -% set the plot style parameters - -plot_style - -outdir = '../../../out/Heskestad_Flame_Height/'; -M = csvread([outdir,'box_height.csv'],1,0); -Qs = M(:,1); -Q = M(:,2); -RI = {'_RI=05','_RI=10','_RI=20'}; -QI = {'p1','p2','p5','1','2','5','10','20','50','100','200','500','1000','2000','5000','10000'}; - -K(4)=loglog(Qs,Q,'k-'); hold on - -for j=1:length(Qs) - for i=1:3 - filename = ['Qs=',QI{j},RI{i},'_hrr.csv']; - A = csvread([outdir,filename],2,0); - n = length(A(:,1)); - Q_fds = mean(A(round(n/2):n,2)); - if (i==1); K(1)=loglog(Qs(j),Q_fds,'ksq'); end - if (i==2); K(2)=loglog(Qs(j),Q_fds,'r^'); end - if (i==3); K(3)=loglog(Qs(j),Q_fds,'go'); end - end -end - -axis([.05 10^4 10^2 10^8]) - -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Label_Font_Size) -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - -xlabel('{\it Q*}','Interpreter',Font_Interpreter,'FontSize',Label_Font_Size) -ylabel('Heat Release Rate (kW)','Interpreter',Font_Interpreter,'FontSize',Label_Font_Size) -text(.1,2e7,'Flame Height Heat Release Verification','FontSize',Title_Font_Size,'FontName',Font_Name,'Interpreter',Font_Interpreter) -hh=legend(K,'{\itD}*/{\it\deltax}=5','{\itD}*/{\it\deltax}=10','{\itD}*/{\it\deltax}=20','correct','Location','Southeast'); -set(hh,'Interpreter',Font_Interpreter,'FontSize',Key_Font_Size) - -% add version if file is available - -Git_Filename = [outdir,'Qs=1_RI=05_git.txt']; -addverstr(gca,Git_Filename,'loglog') - -% print to pdf - -plotdir = '../../Manuals/FDS_Validation_Guide/SCRIPT_FIGURES/Heskestad/'; -Plot_Filename = 'Flame_Height_check_hrr'; - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); -display(['Printing plot Flame_Height_check_hrr.pdf ...']) -print(gcf,'-dpdf',[plotdir,Plot_Filename]) - - - diff --git a/Utilities/Matlab/scripts/compression_wave.m b/Utilities/Matlab/scripts/compression_wave.m deleted file mode 100644 index 8c38cdfebd2..00000000000 --- a/Utilities/Matlab/scripts/compression_wave.m +++ /dev/null @@ -1,198 +0,0 @@ -% McDermott -% 4-8-10 -% compression_wave.m - -close all -clear all - -data_dir = '../../Verification/Scalar_Analytical_Solution/'; -plot_dir = '../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/'; - -skip_case = 0; - -if ~exist([data_dir,'compression_wave_FL0_16_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_16_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL0_32_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_32_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL0_64_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_64_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL0_128_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_128_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end - -if ~exist([data_dir,'compression_wave_FL2_16_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_16_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL2_32_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_32_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL2_64_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_64_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL2_128_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_128_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end - -if ~exist([data_dir,'compression_wave_FL4_16_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_16_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL4_32_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_32_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL4_64_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_64_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end -if ~exist([data_dir,'compression_wave_FL4_128_devc.csv']) - display(['Error: File ' [data_dir,'compression_wave_FL0_128_devc.csv'] ' does not exist. Skipping case.']) - skip_case = 1; -end - -if skip_case - return -end - -% central differencing, FL=0 - -M_FL0_16 = csvread([data_dir,'compression_wave_FL0_16_devc.csv'],2,0); -M_FL0_32 = csvread([data_dir,'compression_wave_FL0_32_devc.csv'],2,0); -M_FL0_64 = csvread([data_dir,'compression_wave_FL0_64_devc.csv'],2,0); -M_FL0_128 = csvread([data_dir,'compression_wave_FL0_128_devc.csv'],2,0); -t_FL0_16 = M_FL0_16(:,1); rho_fds_FL0_16 = M_FL0_16(:,3); -t_FL0_32 = M_FL0_32(:,1); rho_fds_FL0_32 = M_FL0_32(:,3); -t_FL0_64 = M_FL0_64(:,1); rho_fds_FL0_64 = M_FL0_64(:,3); -t_FL0_128 = M_FL0_128(:,1); rho_fds_FL0_128 = M_FL0_128(:,3); - -% Superbee limiter, FL=2 -M_FL2_16 = csvread([data_dir,'compression_wave_FL2_16_devc.csv'],2,0); -M_FL2_32 = csvread([data_dir,'compression_wave_FL2_32_devc.csv'],2,0); -M_FL2_64 = csvread([data_dir,'compression_wave_FL2_64_devc.csv'],2,0); -M_FL2_128 = csvread([data_dir,'compression_wave_FL2_128_devc.csv'],2,0); -t_FL2_16 = M_FL2_16(:,1); rho_fds_FL2_16 = M_FL2_16(:,3); -t_FL2_32 = M_FL2_32(:,1); rho_fds_FL2_32 = M_FL2_32(:,3); -t_FL2_64 = M_FL2_64(:,1); rho_fds_FL2_64 = M_FL2_64(:,3); -t_FL2_128 = M_FL2_128(:,1); rho_fds_FL2_128 = M_FL2_128(:,3); - -% CHARM limiter, FL=4 -M_FL4_16 = csvread([data_dir,'compression_wave_FL4_16_devc.csv'],2,0); -M_FL4_32 = csvread([data_dir,'compression_wave_FL4_32_devc.csv'],2,0); -M_FL4_64 = csvread([data_dir,'compression_wave_FL4_64_devc.csv'],2,0); -M_FL4_128 = csvread([data_dir,'compression_wave_FL4_128_devc.csv'],2,0); -t_FL4_16 = M_FL4_16(:,1); rho_fds_FL4_16 = M_FL4_16(:,3); -t_FL4_32 = M_FL4_32(:,1); rho_fds_FL4_32 = M_FL4_32(:,3); -t_FL4_64 = M_FL4_64(:,1); rho_fds_FL4_64 = M_FL4_64(:,3); -t_FL4_128 = M_FL4_128(:,1); rho_fds_FL4_128 = M_FL4_128(:,3); - -% analytical solution - -a = 2; -c = 3; - -L = 2*pi; -x = 1.5*pi; -y = 1.5*pi; - -rho_FL0_16 = compression_wave_soln(rho_fds_FL0_16(1),x-L/32,y-L/32,a,c,t_FL0_16); error_FL0_16 = norm(rho_fds_FL0_16-rho_FL0_16)/length(t_FL0_16); -rho_FL0_32 = compression_wave_soln(rho_fds_FL0_32(1),x-L/64,y-L/64,a,c,t_FL0_32); error_FL0_32 = norm(rho_fds_FL0_32-rho_FL0_32)/length(t_FL0_32); -rho_FL0_64 = compression_wave_soln(rho_fds_FL0_64(1),x-L/128,y-L/128,a,c,t_FL0_64); error_FL0_64 = norm(rho_fds_FL0_64-rho_FL0_64)/length(t_FL0_64); -rho_FL0_128 = compression_wave_soln(rho_fds_FL0_128(1),x-L/256,y-L/256,a,c,t_FL0_128); error_FL0_128 = norm(rho_fds_FL0_128-rho_FL0_128)/length(t_FL0_128); - -rho_FL2_16 = compression_wave_soln(rho_fds_FL2_16(1),x-L/32,y-L/32,a,c,t_FL2_16); error_FL2_16 = norm(rho_fds_FL2_16-rho_FL2_16)/length(t_FL2_16); -rho_FL2_32 = compression_wave_soln(rho_fds_FL2_32(1),x-L/64,y-L/64,a,c,t_FL2_32); error_FL2_32 = norm(rho_fds_FL2_32-rho_FL2_32)/length(t_FL2_32); -rho_FL2_64 = compression_wave_soln(rho_fds_FL2_64(1),x-L/128,y-L/128,a,c,t_FL2_64); error_FL2_64 = norm(rho_fds_FL2_64-rho_FL2_64)/length(t_FL2_64); -rho_FL2_128 = compression_wave_soln(rho_fds_FL2_128(1),x-L/256,y-L/256,a,c,t_FL2_128); error_FL2_128 = norm(rho_fds_FL2_128-rho_FL2_128)/length(t_FL2_128); - -rho_FL4_16 = compression_wave_soln(rho_fds_FL4_16(1),x-L/32,y-L/32,a,c,t_FL4_16); error_FL4_16 = norm(rho_fds_FL4_16-rho_FL4_16)/length(t_FL4_16); -rho_FL4_32 = compression_wave_soln(rho_fds_FL4_32(1),x-L/64,y-L/64,a,c,t_FL4_32); error_FL4_32 = norm(rho_fds_FL4_32-rho_FL4_32)/length(t_FL4_32); -rho_FL4_64 = compression_wave_soln(rho_fds_FL4_64(1),x-L/128,y-L/128,a,c,t_FL4_64); error_FL4_64 = norm(rho_fds_FL4_64-rho_FL4_64)/length(t_FL4_64); -rho_FL4_128 = compression_wave_soln(rho_fds_FL4_128(1),x-L/256,y-L/256,a,c,t_FL4_128); error_FL4_128 = norm(rho_fds_FL4_128-rho_FL4_128)/length(t_FL4_128); - -figure -plot_style -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - -H(5)=plot(t_FL4_128,rho_FL4_128,'k-','LineWidth',Line_Width); hold on -H(1)=plot(t_FL4_16,rho_fds_FL4_16,'c--','LineWidth',Line_Width); hold on -H(2)=plot(t_FL4_32,rho_fds_FL4_32,'g--','LineWidth',Line_Width); hold on -H(3)=plot(t_FL4_64,rho_fds_FL4_64,'b--','LineWidth',Line_Width); hold on -H(4)=plot(t_FL4_128,rho_fds_FL4_128,'r--','LineWidth',Line_Width); hold on - -xlabel('Time (s)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter,'FontName',Font_Name) -ylabel('Density (kg/m³)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter,'FontName',Font_Name) -axis([0 12.5 0 8]) -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Title_Font_Size) -legend_handle=legend(H,'FDS N=16','FDS N=32','FDS N=64','FDS N=128','Analytical Solution','Location','NorthEast'); -set(legend_handle,'Interpreter',Font_Interpreter); -set(legend_handle,'Fontsize',Key_Font_Size); -set(legend_handle,'Box','on'); - -% add Git revision if file is available - -Git_Filename = [data_dir,'compression_wave_FL0_16_git.txt']; -addverstr(gca,Git_Filename,'linear') - -% print to pdf -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); -print(gcf,'-dpdf',[plot_dir,'compression_wave_time_series']) - -% convergence plot - -figure -plot_style -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - -h = 2*pi./[16 32 64 128]; -e_FL0 = [error_FL0_16 error_FL0_32 error_FL0_64 error_FL0_128]; -e_FL2 = [error_FL2_16 error_FL2_32 error_FL2_64 error_FL2_128]; -e_FL4 = [error_FL4_16 error_FL4_32 error_FL4_64 error_FL4_128]; -H(1)=loglog(h,e_FL0,'b*-','LineWidth',Line_Width); hold on -H(2)=loglog(h,e_FL2,'ro-','LineWidth',Line_Width); hold on -H(3)=loglog(h,e_FL4,'g^-','LineWidth',Line_Width); hold on -H(4)=loglog(h,.1*h,'k--','LineWidth',Line_Width); -H(5)=loglog(h,.1*h.^2,'k-','LineWidth',Line_Width); - -xlabel('Grid Spacing (m)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter,'FontName',Font_Name) -ylabel('L2 Error (kg/m³)','FontSize',Title_Font_Size,'Interpreter',Font_Interpreter,'FontName',Font_Name) -axis([1e-2 1e0 1e-4 1e-1]) -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Title_Font_Size) -legend_handle=legend(H(1:5),'Central','Superbee','CHARM','{\it O}({\it\deltax})','{\it O}({\it\deltax^2})','Location','NorthWest'); -set(legend_handle,'Interpreter',Font_Interpreter); -set(legend_handle,'Fontsize',Key_Font_Size); -set(legend_handle,'Box','on'); - -% add Git revision if file is available - -Git_Filename = [data_dir,'compression_wave_FL0_16_git.txt']; -addverstr(gca,Git_Filename,'loglog') - -% print to pdf -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); -print(gcf,'-dpdf',[plot_dir,'compression_wave_convergence']) - - - diff --git a/Utilities/Matlab/scripts/compression_wave_soln.m b/Utilities/Matlab/scripts/compression_wave_soln.m deleted file mode 100644 index ec4713d0ae8..00000000000 --- a/Utilities/Matlab/scripts/compression_wave_soln.m +++ /dev/null @@ -1,22 +0,0 @@ -% McDermott -% 4-8-10 -% compression_wave_soln.m - -function [rho] = compression_wave_soln(rho0,x,y,a,c,t); - -b = sqrt(-1+a^2); -d = sqrt(-1+c^2); - -x0 = 2*atan( b/a*tan( atan( (1+a*tan(x/2))/b ) - b*t/2 ) - 1/a ); -y0 = 2*atan( d/c*tan( atan( (1+c*tan(y/2))/d ) - d*t/2 ) - 1/c ); - -Ix0 = log( -a^2 - cos( 2*atan( (1+a*tan(x0/2))/b ) ) + b*sin( 2*atan( (1+a*tan(x0/2))/b ) ) ); -Iy0 = log( -c^2 - cos( 2*atan( (1+c*tan(y0/2))/d ) ) + d*sin( 2*atan( (1+c*tan(y0/2))/d ) ) ); - -Ix = log( -a^2 - cos(b*t + 2*atan( (1+a*tan(x0/2))/b ) ) + b*sin(b*t + 2*atan( (1+a*tan(x0/2))/b ) ) ); -Iy = log( -c^2 - cos(d*t + 2*atan( (1+c*tan(y0/2))/d ) ) + d*sin(d*t + 2*atan( (1+c*tan(y0/2))/d ) ) ); - -q0 = log(rho0); -q = q0 + Ix-Ix0 + Iy-Iy0; - -rho = exp(q); \ No newline at end of file diff --git a/Utilities/Matlab/scripts/flat_fire_comparison.m b/Utilities/Matlab/scripts/flat_fire_comparison.m deleted file mode 100644 index 710261ebfce..00000000000 --- a/Utilities/Matlab/scripts/flat_fire_comparison.m +++ /dev/null @@ -1,39 +0,0 @@ -% Trettel -% flat fire comparison -% This script just computes the exact solution. The actual plots are made by dataplot, and the exact solution is committed to the repo. - -close all -clear all - -V_0 = 400; -h = 8; -g = 9.8; -C_d = 0.2; -rho_a = 1.2; -rho_d = 1000; -D = 5e-3; - -K = 3 * rho_a * C_d / (4 * rho_d * D); - -dt = 0.05; -tend = 1.65; - -tvec = 0:dt:tend; - -xexact = log(V_0 * K * tvec + 1) / K; -yexact = h... - + (g / (2 * (K * V_0)^2)) * log(V_0 * K * tvec + 1)... - - g * tvec .^ 2 / 4 ... - - g * tvec ./ (2 * K * V_0); -uexact = V_0 ./ (V_0 * K * tvec + 1); -vexact = g ./ (2 * K * V_0 * (K * V_0 * tvec + 1))... - - g * tvec ./ 2 ... - - g / (2 * K * V_0); - -mat(:,1) = tvec; -mat(:,2) = xexact -mat(:,3) = yexact -mat(:,4) = uexact -mat(:,5) = vexact -csvwrite('flat_fire.csv',mat) - diff --git a/Utilities/Matlab/scripts/geom_channel_test.m b/Utilities/Matlab/scripts/geom_channel_test.m deleted file mode 100644 index 736164a9181..00000000000 --- a/Utilities/Matlab/scripts/geom_channel_test.m +++ /dev/null @@ -1,27 +0,0 @@ -% McDermott -% 8-11-2021 -% geom_channel_test.m -% -% Temporary home for checking AREA integration for geom_channel.fds - -close all -clear all - -M = importdata('../../Verification/Complex_Geometry/geom_channel_devc.csv',',',2); - -U0 = M.data(end,find(strcmp(M.colheaders,'"U0"'))); -A0 = M.data(end,find(strcmp(M.colheaders,'"A0"'))); -U1 = M.data(end,find(strcmp(M.colheaders,'"U1"'))); -A1 = M.data(end,find(strcmp(M.colheaders,'"A1"'))); -U2 = M.data(end,find(strcmp(M.colheaders,'"U2"'))); -A2 = M.data(end,find(strcmp(M.colheaders,'"A2"'))); -U3 = M.data(end,find(strcmp(M.colheaders,'"U3"'))); -A3 = M.data(end,find(strcmp(M.colheaders,'"A3"'))); - -% check errors -error_tolerance = 1e-6; - -A0_error = abs(A0-1); if A0_error>error_tolerance; display(['Matlab Warning: geom_channel.fds A0_error = ',num2str(A0_error)]); end -A1_error = abs(A1-1); if A1_error>error_tolerance; display(['Matlab Warning: geom_channel.fds A1_error = ',num2str(A1_error)]); end -A2_error = abs(A2-1); if A2_error>error_tolerance; display(['Matlab Warning: geom_channel.fds A2_error = ',num2str(A2_error)]); end -A3_error = abs(A3-1); if A3_error>error_tolerance; display(['Matlab Warning: geom_channel.fds A3_error = ',num2str(A3_error)]); end \ No newline at end of file diff --git a/Utilities/Matlab/scripts/geom_positive_errors.m b/Utilities/Matlab/scripts/geom_positive_errors.m deleted file mode 100644 index 8327227e363..00000000000 --- a/Utilities/Matlab/scripts/geom_positive_errors.m +++ /dev/null @@ -1,51 +0,0 @@ -% Emanuele Gissi -% 6-7-2017 -% geom_positive_errors.m - -% This script catches errors produced during FDS setup of bad geometries cases. -% The success condition is the existance of the error message. The error message has -% its 'ERROR' label replaced by 'SUCCESS' label, thanks to the MISC parameter -% POSITIVE_ERROR_TEST=.TRUE. in the bad geometries cases. -% The absence of the 'ERROR' label makes current Firebot to consider the run a success. - -% Clean up - -close all -clear all - -% Cases and error strings - -infile{1} = '../../Verification/Complex_Geometry/geom_bad_inconsistent_normals.err'; -errstring{1} = 'Face normals are probably pointing in the wrong direction'; - -infile{2} = '../../Verification/Complex_Geometry/geom_bad_open_surface.err'; -errstring{2} = 'Open geometry at edge'; - -infile{3} = '../../Verification/Complex_Geometry/geom_bad_non_manifold_edge.err'; -errstring{3} = 'Non manifold geometry in adjacent faces at edge'; - -infile{4} = '../../Verification/Complex_Geometry/geom_bad_inverted_normals.err'; -errstring{4} = 'Face normals are probably pointing in the wrong direction'; - -% FIXME This error condition is currently not caught -% It will be caught when boolean operations are performed -%infile{7} = '../../Verification/Complex_Geometry/geom_bad_non_manifold_vert.err'; -%errstring{7} = 'Non manifold geometry at vertex' - -% Check existance of errstring in each infile - -for n = 1:4 - % infile exists? - if ~exist(infile{n},'file') - display(['Error: File ',infile{n},' does not exist. Skipping case.']) - continue - end - - % errstring in infile? - errfile = fileread(infile{n}); - if isempty(strfind(errfile, errstring{n})); - display(['Error: File ',infile{n},':']) - display([' Does not contain the following positive error message:']) - display([' <',errstring{n},'>']) - end -end diff --git a/Utilities/Matlab/scripts/mass_balance_gas_volume.m b/Utilities/Matlab/scripts/mass_balance_gas_volume.m deleted file mode 100644 index 2dd0127638e..00000000000 --- a/Utilities/Matlab/scripts/mass_balance_gas_volume.m +++ /dev/null @@ -1,71 +0,0 @@ -% McDermott -% 2 Dec 2021 -% mass_balance_gas_volume.m - -close all -clear all - -plot_mass_balance('mass_balance_gas_volume',''); - -function [] = plot_mass_balance(chid,title_text) - -plot_style -figure -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - -ddir = '../../Verification/Species/'; -M = importdata([ddir,chid,'_devc.csv'],',',2); - -t = M.data(:,1); -m = M.data(:,find(strcmp(M.colheaders,'M'))); -dmdt = zeros(length(t),1); -for i=2:length(t) - dmdt(i) = (m(i)-m(i-1))/(t(i)-t(i-1)); -end - -mf_x1 = M.data(:,find(strcmp(M.colheaders,'TMF_X1'))); -mf_x2 = M.data(:,find(strcmp(M.colheaders,'TMF_X2'))); -mf_y1 = M.data(:,find(strcmp(M.colheaders,'TMF_Y1'))); -mf_y2 = M.data(:,find(strcmp(M.colheaders,'TMF_Y2'))); -mf_z1 = M.data(:,find(strcmp(M.colheaders,'TMF_Z1'))); -mf_z2 = M.data(:,find(strcmp(M.colheaders,'TMF_Z2'))); - -bal = dmdt + (mf_x2-mf_x1 + mf_y2-mf_y1 + mf_z2-mf_z1); - -plot(t,zeros(1,length(t)),'k-'); hold on -H(1)=plot(t,dmdt); -H(2)=plot(t,(mf_x2-mf_x1 + mf_y2-mf_y1 + mf_z2-mf_z1)); -H(3)=plot(t,bal); - -ylabel('mass flow (kg/s)', 'FontSize',Label_Font_Size) -xlabel('time (s)', 'FontSize',Label_Font_Size) - -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Label_Font_Size) - -lh=legend(H,'dm/dt','out-in','dm/dt+out-in'); -set(lh,'FontSize',Key_Font_Size) - -text(100,18e-3,title_text,'FontSize',Title_Font_Size,'FontName',Font_Name,'Interpreter',Font_Interpreter) - -% check balance and report error - -mass_error = max(abs(bal)); -if mass_error > 1e-6 - disp(['Matlab Warning: mass error = ',num2str(mass_error),' in ',chid]) -end - -% add version string if file is available - -Git_Filename = [ddir,chid,'_git.txt']; -addverstr(gca,Git_Filename,'linear') - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); -print(gcf,'-dpdf',['../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/',chid]) - -end diff --git a/Utilities/Matlab/scripts/mass_balance_reac.m b/Utilities/Matlab/scripts/mass_balance_reac.m deleted file mode 100644 index e8b4f2ef2ba..00000000000 --- a/Utilities/Matlab/scripts/mass_balance_reac.m +++ /dev/null @@ -1,83 +0,0 @@ -% McDermott -% 10 July 2020 -% mass_balance_reac.m - -close all -clear all - -plot_mass_balance('mass_balance_reac','Propane Mass Balance','PROPANE','C3H8'); -plot_mass_balance('mass_balance_reac','Oxygen Mass Balance','OXYGEN','O2'); -plot_mass_balance('mass_balance_reac','Nitrogen Mass Balance','NITROGEN','N2'); -plot_mass_balance('mass_balance_reac','Carbon Dioxide Mass Balance','CARBON DIOXIDE','CO2'); -plot_mass_balance('mass_balance_reac','Water Vapor Mass Balance','WATER VAPOR','H2O'); -plot_mass_balance('mass_balance_reac','Soot Mass Balance','SOOT','Soot'); - -function [] = plot_mass_balance(chid,title_text,mass_id,devc_id) - -plot_style -figure -set(gca,'Units',Plot_Units) -set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - -ddir = '../../Verification/Species/'; -M = importdata([ddir,chid,'_mass.csv'],',',2); - -t = M.data(:,1); -m = M.data(:,find(strcmp(M.colheaders,mass_id))); -dmdt = zeros(length(t),1); -for i=2:length(t) - dmdt(i) = (m(i)-m(i-1))/(t(i)-t(i-1)); -end - -F = importdata([ddir,chid,'_devc.csv'],',',2); - -mdot_out = F.data(:,find(strcmp(F.colheaders,[devc_id,' xmax']))) ... - + F.data(:,find(strcmp(F.colheaders,[devc_id,' xmin']))) ... - + F.data(:,find(strcmp(F.colheaders,[devc_id,' ymin']))) ... - + F.data(:,find(strcmp(F.colheaders,[devc_id,' ymax']))) ... - + F.data(:,find(strcmp(F.colheaders,[devc_id,' Burner']))); - -gen = F.data(:,find(strcmp(F.colheaders,[devc_id,' mdot reac']))); - -bal = -dmdt + mdot_out + gen; - -%plot(t,zeros(1,length(t)),'k-'); hold on -H(1)=plot(t,dmdt,'g-'); hold on -H(2)=plot(t,mdot_out,'b-'); -H(3)=plot(t,gen,'r-'); -H(4)=plot(t,bal,'k-'); - -ylabel('Mass Flow (kg/s)', 'FontSize',Label_Font_Size) -xlabel('Time (s)', 'FontSize',Label_Font_Size) - -set(gca,'FontName',Font_Name) -set(gca,'FontSize',Label_Font_Size) - -lh=legend(H,'accumulation','in - out','generation', 'balance','location','southwest'); -set(lh,'FontSize',Key_Font_Size) - -xl=xlim; -yl=ylim; - -text(xl(1)+0.05*(xl(2)-xl(1)),yl(1)+0.9*(yl(2)-yl(1)),title_text,'FontSize',Title_Font_Size,'FontName',Font_Name,'Interpreter',Font_Interpreter) - -% check balance and report error - -mass_error = abs(mean(bal))/m(end); % error normalized by the total mass in the domain -if mass_error > 1.5e-3 - disp(['Matlab Warning: mass error = ',num2str(mass_error),' in ',chid,' for species ',mass_id]) -end - -% add version string if file is available - -Git_Filename = [ddir,chid,'_git.txt']; -addverstr(gca,Git_Filename,'linear') - -set(gcf,'Visible',Figure_Visibility); -set(gcf,'Units',Paper_Units); -set(gcf,'PaperUnits',Paper_Units); -set(gcf,'PaperSize',[Paper_Width Paper_Height]); -set(gcf,'Position',[0 0 Paper_Width Paper_Height]); -print(gcf,'-dpdf',['../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/',chid,'_',devc_id]) - -end diff --git a/Utilities/Matlab/scripts/reaction_rate.m b/Utilities/Matlab/scripts/reaction_rate.m deleted file mode 100644 index 15428f60047..00000000000 --- a/Utilities/Matlab/scripts/reaction_rate.m +++ /dev/null @@ -1,25 +0,0 @@ -% McGrattan -% 8-10-09 -% reaction_rate.m -% -% input: -% T = temperature (K) -% Y = mass fraction -% -% output: -% r = reaction rate - -function r = reaction_rate(T,Y) - -global dTdt -global R0 -global E -global A -global residue - -r(1) = -A(1).*Y(1)*exp(-E(1)./(R0.*T))./dTdt; -r(2) = -A(2).*Y(2)*exp(-E(2)./(R0.*T))./dTdt; -r(3) = -residue(2)*r(2); -r=r'; - - diff --git a/Utilities/Matlab/scripts/rotcube_cc_mms_error.m b/Utilities/Matlab/scripts/rotcube_cc_mms_error.m deleted file mode 100644 index cb19fd23ed5..00000000000 --- a/Utilities/Matlab/scripts/rotcube_cc_mms_error.m +++ /dev/null @@ -1,345 +0,0 @@ -%!/usr/bin/matlab -%Vanella rotcube_cc_mms_error.m: -%05-24-2019 - -close all -clear all - -% Rotated cube parameters: -mu = 0.01; -D = 0.01; -rho= 1.00; -L = pi; % Cube side length. - -nwave=1.; % Wave number on analytical solution. -gam =pi/2; -Az =0.1; -meanz=0.15; -displ =pi; -dispxy=-1/2*displ*ones(2,1); - -% Analytical solutions: -ann_z = @(x,y,t) ... -Az/3.*sin(t)*(1.-cos(2.*nwave*(x-gam)))*(1.-cos(2.*nwave*(y-gam))) - ... -Az/3.*sin(t)+meanz; - -ann_u = @(x,y,t) ... --sin(t)*sin(nwave*x)^2 * sin(2.*nwave*y); - -ann_v = @(x,y,t) ... -sin(t)*sin(2.*nwave*x) * sin(nwave*y)^2; - -ann_p = @(x,y,t) ... -sin(t)/4.*(2.+cos(2.*nwave*x))*(2.+cos(2.*nwave*y))-sin(t); - -ann_H = @(x,y,t) ... -(ann_u(x,y,t)^2+ann_v(x,y,t)^2)/2 + ann_p(x,y,t)/rho; - -datadir = '../../Verification/Complex_Geometry/'; - -ifile_s=1; -ifile_f=3; -ifile_ang = [0; atan(1./2.); atan(1.)]; -ifile_str = {'0deg_';'27deg_';'45deg_'}; -jfile_str = {'stm';'obs'}; - -skip_case = 0; -nfile = 0; -for ifile=ifile_s:ifile_f - - jfile_s = 1; - jfile_f = 1; - if(ifile==1) - jfile_f = 2; - end - for jfile=jfile_s:jfile_f - nfile = nfile+1; - file(nfile).name = { ['rotated_cube_' ifile_str{ifile} '32_' jfile_str{jfile} '_mms.csv'], ... - ['rotated_cube_' ifile_str{ifile} '64_' jfile_str{jfile} '_mms.csv'], ... - ['rotated_cube_' ifile_str{ifile} '128_' jfile_str{jfile} '_mms.csv'], ... - ['rotated_cube_' ifile_str{ifile} '256_' jfile_str{jfile} '_mms.csv'], ...}; - ['rotated_cube_' ifile_str{ifile} '320_' jfile_str{jfile} '_mms.csv']}; - file(nfile).nameout = ['rotated_cube_' ifile_str{ifile} jfile_str{jfile} '_mms_convergence']; - file(nfile).rotang=ifile_ang(ifile); - for n=1:length(file(nfile).name) - if ~exist([datadir,file(nfile).name{n}]) - display(['Error: File ' [datadir,file(nfile).name{n}] ... - ' does not exist. Skipping case.']) - skip_case = 1; - end - end - end -end -if skip_case - return -end - -% Compute errors from *_mms.csv files: -for ifile=1:nfile - - e_z = zeros(length(file(ifile).name),1); - e_u = zeros(length(file(ifile).name),1); - e_v = zeros(length(file(ifile).name),1); - e_H = zeros(length(file(ifile).name),1); - - rotangle = file(ifile).rotang; - ROTMAT = [cos(rotangle) -sin(rotangle); sin(rotangle) cos(rotangle)]; - TROTMAT = ROTMAT'; - - for n=1:length(file(ifile).name) - -% disp(file(ifile).name{n}) - - M = importdata([datadir,file(ifile).name{n}],',',1); % FDS results. - vec = str2num(M.textdata{1,1}); % exact time value - - ntot_u = vec(1); - ntot_v = vec(2); - ntot_c = vec(3); - T = vec(4); - dx = vec(5); - dy = vec(6); - - file(ifile).dx(n) = dx; - - % inialize error arrays - e_z_vec = zeros(ntot_c,1); - e_u_vec = zeros(ntot_u,1); - e_v_vec = zeros(ntot_v,1); - e_H_vec = zeros(ntot_c,1); - - - p = 0; - % First u: - Atot_u = 0; - for i=1:ntot_u - p=p+1; - iscf = M.data(p,1); - xglob= M.data(p,2); - yglob= M.data(p,3); - area = M.data(p,4); - uglob= M.data(p,5); - - % Compute analytical value of u in global coordinate system: - XGLOB = [xglob yglob]' - displ; - XLOC = TROTMAT*XGLOB - dispxy; - - % Analytical Velocity in local axes: - ULOC(1,1) = ann_u(XLOC(1),XLOC(2),T); - ULOC(2,1) = ann_v(XLOC(1),XLOC(2),T); - - % Analytical Velocity in global axes: - UGLOB = ROTMAT*ULOC; - - uglob_ann = UGLOB(1); - - % Error: - e_u_vec(i) = uglob - uglob_ann; - - Atot_u = Atot_u + area; - end - - file(ifile).e_u_1(n) = norm(e_u_vec(1:ntot_u),1)/ntot_u; - file(ifile).e_u_2(n) = norm(e_u_vec(1:ntot_u),2)/sqrt(ntot_u); - file(ifile).e_u_i(n) = norm(e_u_vec(1:ntot_u),inf); - - % Second v: - Atot_v = 0; - for i=1:ntot_v - p=p+1; - iscf = M.data(p,1); - xglob= M.data(p,2); - yglob= M.data(p,3); - area = M.data(p,4); - vglob= M.data(p,5); - - % Compute analytical value of v in global coordinate system: - XGLOB = [xglob yglob]' - displ; - XLOC = TROTMAT*XGLOB - dispxy; - - % Analytical Velocity in local axes: - ULOC(1,1) = ann_u(XLOC(1),XLOC(2),T); - ULOC(2,1) = ann_v(XLOC(1),XLOC(2),T); - - % Analytical Velocity in global axes: - UGLOB = ROTMAT*ULOC; - - vglob_ann = UGLOB(2); - - % Error: - e_v_vec(i) = vglob - vglob_ann; - - Atot_v = Atot_v + area; - end - - file(ifile).e_v_1(n) = norm(e_v_vec(1:ntot_v),1)/ntot_v; - file(ifile).e_v_2(n) = norm(e_v_vec(1:ntot_v),2)/sqrt(ntot_v); - file(ifile).e_v_i(n) = norm(e_v_vec(1:ntot_v),inf); - - % Finally cell centered variables: - Vtot_c = 0; - DELP = 0; - p_c = p; - for i=1:ntot_c - p=p+1; - iscf = M.data(p,1); - xglob= M.data(p,2); - yglob= M.data(p,3); - vol = M.data(p,4); - z = M.data(p,5); - H = M.data(p,6); - pres = M.data(p,7); - - % Compute analytical scalars in global coordinate system: - XGLOB = [xglob yglob]' - displ; - XLOC = TROTMAT*XGLOB - dispxy; - - % Analytical Velocity in local axes: - z_ann = ann_z(XLOC(1),XLOC(2),T); - H_ann = ann_H(XLOC(1),XLOC(2),T); - p_ann = ann_p(XLOC(1),XLOC(2),T); - - if(i==1) - DELP = pres - p_ann; - end - - % Error: - e_z_vec(i) = z - z_ann; - e_H_vec(i) = H - H_ann - DELP/rho; - e_p_vec(i) = pres - p_ann - DELP; - - Vtot_c = Vtot_c + vol; - end - - file(ifile).e_z_1(n) = norm(e_z_vec(1:ntot_c),1)/ntot_c; - file(ifile).e_z_2(n) = norm(e_z_vec(1:ntot_c),2)/sqrt(ntot_c); - file(ifile).e_z_i(n) = norm(e_z_vec(1:ntot_c),inf); - - file(ifile).e_H_1(n) = norm(e_H_vec(1:ntot_c),1)/ntot_c; - file(ifile).e_H_2(n) = norm(e_H_vec(1:ntot_c),2)/sqrt(ntot_c); - file(ifile).e_H_i(n) = norm(e_H_vec(1:ntot_c),inf); - - file(ifile).e_p_1(n) = norm(e_p_vec(1:ntot_c),1)/ntot_c; - file(ifile).e_p_2(n) = norm(e_p_vec(1:ntot_c),2)/sqrt(ntot_c); - file(ifile).e_p_i(n) = norm(e_p_vec(1:ntot_c),inf); - [dummy,i] = max(abs(e_p_vec(1:ntot_c))); - file(ifile).e_p_i_loc(n) = i; - file(ifile).e_p_i_xy(n,1:3) = [ M.data(p,7) M.data(p_c+i,2) M.data(p_c+i,3) ]; - - end - - % Error orders: -% p_L2_u = log(file(ifile).e_u_2(1:end-1)./file(ifile).e_u_2(2:end))./... -% log(file(ifile).dx(1:end-1) ./file(ifile).dx(2:end)); -% p_Li_u = log(file(ifile).e_u_i(1:end-1)./file(ifile).e_u_i(2:end))./... -% log(file(ifile).dx(1:end-1) ./file(ifile).dx(2:end)); -% disp(' ') -% disp(['L2 Err u p:' file(ifile).nameout]) -% disp(' dx L2e p2 Linfe pinf') -% for ip=1:length(p_L2_u) -% disp([num2str(file(ifile).dx(ip+1)) ', ' ... -% num2str(file(ifile).e_u_2(ip+1)) ', ' ... -% num2str(p_L2_u(ip)) ', ' ... -% num2str(file(ifile).e_u_i(ip+1)) ', ' ... -% num2str(p_Li_u(ip))]) -% end -% disp(' ') -% p_L2_z = log(file(ifile).e_z_2(1:end-1)./file(ifile).e_z_2(2:end))./... -% log(file(ifile).dx(1:end-1) ./file(ifile).dx(2:end)); -% p_Li_z = log(file(ifile).e_z_i(1:end-1)./file(ifile).e_z_i(2:end))./... -% log(file(ifile).dx(1:end-1) ./file(ifile).dx(2:end)); -% disp(' ') -% disp(['L2 Err z p:' file(ifile).nameout]) -% disp(' dx L2e p2 Linfe pinf') -% for ip=1:length(p_L2_z) -% disp([num2str(file(ifile).dx(ip+1)) ', ' ... -% num2str(file(ifile).e_z_2(ip+1)) ', ' ... -% num2str(p_L2_z(ip)) ', ' ... -% num2str(file(ifile).e_z_i(ip+1)) ', ' ... -% num2str(p_Li_z(ip))]) -% end -% disp(' ') - - plot_style - figure - set(gca,'Units',Plot_Units) - set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height]) - - hh(1)=loglog(file(ifile).dx,file(ifile).e_z_2,'ks-'); - hold on - hh(2)=loglog(file(ifile).dx,file(ifile).e_u_2,'k>-'); - hh(3)=loglog(file(ifile).dx,file(ifile).e_H_2,'k+-'); - hh(4)=loglog(file(ifile).dx,10^-1*file(ifile).dx,'k--'); - hh(5)=loglog(file(ifile).dx,2*10^-3*file(ifile).dx.^2,'k-'); - if(ifile <= 2) - axis([2*10^-3 3*10^-1 10^-7 10^-1]) - elseif(ifile <=3) - axis([10^-2 5*10^-1 10^-7 7*10^-1]) - else - axis([5*10^-3 5*10^-1 10^-7 10^-1]) - end - set(gca,'FontName',Font_Name) - set(gca,'FontSize',Title_Font_Size) - - xlabel('{\it \Deltax} (m)','FontSize',Title_Font_Size,'Interpreter',... - Font_Interpreter,'Fontname',Font_Name) - ylabel('{\it L_2} Error','FontSize',Title_Font_Size,'Interpreter', ... - Font_Interpreter,'Fontname',Font_Name) - lh=legend(hh,'FDS {\it Z}','FDS {\it u}',... - 'FDS {\it H}','{\it O(\Deltax)}','{\it O(\Deltax^2)}',... - 'location','northwest'); - set(lh,'FontSize',Title_Font_Size,'Interpreter',Font_Interpreter,... - 'Fontname',Font_Name) - - % add Git version if file is available - Git_Filename = [datadir,'rotated_cube_45deg_256_stm_git.txt']; - addverstr(gca,Git_Filename,'loglog') - - % print to pdf - set(gcf,'Visible',Figure_Visibility); - set(gcf,'Units',Paper_Units); - set(gcf,'PaperUnits',Paper_Units); - set(gcf,'PaperSize',[Paper_Width Paper_Height]); - set(gcf,'Position',[0 0 Paper_Width Paper_Height]); - - strng=['../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/' ... - file(ifile).nameout]; - - print(gcf,'-dpdf',strng) - - if (ifile==2) % OBST case: - % check errors - if file(ifile).e_z_2(end) > 2e-6 - display(['Matlab Warning: Neumann BC Species in rotated_cube OBST ' ... - file(ifile).name{end} ' is out of tolerance. e_z = ', ... - num2str(file(ifile).e_z_2(end))]) - end - if file(ifile).e_u_2(end) > 2e-5 - display(['Matlab Warning: Velocity in rotated_cube OBST ' ... - file(ifile).name{end} ' is out of tolerance. e_u = ', ... - num2str(file(ifile).e_u_2(end))]) - end - if file(ifile).e_H_2(end) > 1.5e-3 - display(['Matlab Warning: Pressure in rotated_cube OBST ' ... - file(ifile).name{end} ' is out of tolerance. e_H = ', ... - num2str(file(ifile).e_H_2(end))]) - end - elseif(ifile==3) - % check errors - if file(ifile).e_z_2(end) > 6.5e-5 - display(['Matlab Warning: Neumann BC Species in rotated_cube 27deg stm ' ... - file(ifile).name{end} ' is out of tolerance. e_z = ', ... - num2str(file(ifile).e_z_2(end))]) - end - if file(ifile).e_u_2(end) > 8.5e-4 - display(['Matlab Warning: Velocity in rotated_cube 27deg stm ' ... - file(ifile).name{end} ' is out of tolerance. e_u = ', ... - num2str(file(ifile).e_u_2(end))]) - end - if file(ifile).e_H_2(end) > 2.5e-3 - display(['Matlab Warning: Pressure in rotated_cube 27deg stm ' ... - file(ifile).name{end} ' is out of tolerance. e_H = ', ... - num2str(file(ifile).e_H_2(end))]) - end - end -end diff --git a/Utilities/Python/FDS_validation_script.py b/Utilities/Python/FDS_validation_script.py index 15a5d8dc0d6..24981b6f618 100644 --- a/Utilities/Python/FDS_validation_script.py +++ b/Utilities/Python/FDS_validation_script.py @@ -7,14 +7,12 @@ # Scripts to run prior to dataplot -print("catchpole_spread_rates..."); runpy.run_path("./scripts/catchpole_spread_rates.py", run_name="__main__") print("NIST_deposition_gauge..."); runpy.run_path("./scripts/NIST_deposition_gauge.py", run_name="__main__") print("flame_height..."); runpy.run_path("./scripts/flame_height.py", run_name="__main__") print("NIST_RSE..."); runpy.run_path("./scripts/NIST_RSE.py", run_name="__main__") print("sippola_aerosol_deposition..."); runpy.run_path("./scripts/sippola_aerosol_deposition.py", run_name="__main__") print("layer_height..."); runpy.run_path("./scripts/layer_height.py", run_name="__main__") print("NIST_NRC_Corner_Effects..."); runpy.run_path("./scripts/NIST_NRC_Corner_Effects.py", run_name="__main__") -# print("fm_data_center..."); runpy.run_path("./scripts/fm_data_center.py", run_name="__main__") print("LNG_Dispersion..."); runpy.run_path("./scripts/LNG_Dispersion.py", run_name="__main__") print("LNG_wind_profiles..."); runpy.run_path("./scripts/LNG_wind_profiles.py", run_name="__main__") print("FM_Vertical_Wall_Flames..."); runpy.run_path("./scripts/FM_Vertical_Wall_Flames.py", run_name="__main__") @@ -61,8 +59,10 @@ # Special cases +print("Backward_Facing_Step..."); runpy.run_path("./scripts/Backward_Facing_Step.py", run_name="__main__") print("Beyler_Hood..."); runpy.run_path("./scripts/Beyler_Hood.py", run_name="__main__") print("BRE_LEMTA_Sprays..."); runpy.run_path("./scripts/BRE_LEMTA_Sprays.py", run_name="__main__") +print("catchpole_spread_rates..."); runpy.run_path("./scripts/catchpole_spread_rates.py", run_name="__main__") print("FHWA_Tunnel..."); runpy.run_path("./scripts/FHWA_Tunnel.py", run_name="__main__") print("FM_FPRF_Datacenter..."); runpy.run_path("./scripts/FM_FPRF_Datacenter.py", run_name="__main__") print("Heskestad_Flame_Height_2..."); runpy.run_path("./scripts/Heskestad_Flame_Height_2.py", run_name="__main__") diff --git a/Utilities/Python/FDS_verification_script.py b/Utilities/Python/FDS_verification_script.py index c3610b7eb47..70fe004acf9 100644 --- a/Utilities/Python/FDS_verification_script.py +++ b/Utilities/Python/FDS_verification_script.py @@ -43,11 +43,14 @@ print("atmospheric_boundary_layer..."); runpy.run_path("./scripts/atmospheric_boundary_layer.py", run_name="__main__") print("blasius..."); runpy.run_path("./scripts/blasius.py", run_name="__main__") +print('compression_wave...'); runpy.run_path("./scripts/compression_wave.py", run_name="__main__") print('extinction...'); runpy.run_path("./scripts/extinction.py", run_name="__main__") print("fan_curve..."); runpy.run_path("./scripts/fan_curve.py", run_name="__main__") print("favre_test..."); runpy.run_path("./scripts/favre_test.py", run_name="__main__") print("fds_moody_chart..."); runpy.run_path("./scripts/fds_moody_chart.py", run_name="__main__") print("fluid_part..."); runpy.run_path("./scripts/fluid_part.py", run_name="__main__") +print("geom_channel_test..."); runpy.run_path("./scripts/geom_channel_test.py", run_name="__main__") +print("geom_positive_errors..."); runpy.run_path("./scripts/geom_positive_errors.py", run_name="__main__") print("heated_channel..."); runpy.run_path("./scripts/heated_channel.py", run_name="__main__") print('hot_layer_collapse...'); runpy.run_path("./scripts/hot_layer_collapse.py", run_name="__main__") print("ht3d_sphere..."); runpy.run_path("./scripts/ht3d_sphere.py", run_name="__main__") @@ -57,6 +60,8 @@ print("level_set_ellipse..."); runpy.run_path("./scripts/level_set_ellipse.py", run_name="__main__") print("mesh_transformation..."); runpy.run_path("./scripts/mesh_transformation.py", run_name="__main__") print("mass_balance..."); runpy.run_path("./scripts/mass_balance.py", run_name="__main__") +print("mass_balance_gas_volume..."); runpy.run_path("./scripts/mass_balance_gas_volume.py", run_name="__main__") +print("mass_balance_reac..."); runpy.run_path("./scripts/mass_balance_reac.py", run_name="__main__") print("openmp_timing_benchmarks..."); runpy.run_path("./scripts/openmp_timing_benchmarks.py", run_name="__main__") print("particle_size_distribution..."); runpy.run_path("./scripts/particle_size_distribution.py", run_name="__main__") print("plate_view_factor..."); runpy.run_path("./scripts/plate_view_factor.py", run_name="__main__") @@ -64,6 +69,7 @@ print("pyrolysis..."); runpy.run_path("./scripts/pyrolysis.py", run_name="__main__") print('radiating_polygon...'); runpy.run_path("./scripts/radiating_polygon.py", run_name="__main__") print("ribbed_channel..."); runpy.run_path("./scripts/ribbed_channel.py", run_name="__main__") +print('rotcube_cc_mms_error...'); runpy.run_path("./scripts/rotcube_cc_mms_error.py", run_name="__main__") print("saad_mms..."); runpy.run_path("./scripts/saad_mms_temporal_error.py", run_name="__main__") print("scaling_tests..."); runpy.run_path("./scripts/scaling_tests.py", run_name="__main__") print("shunn_mms..."); runpy.run_path("./scripts/shunn_mms.py", run_name="__main__") diff --git a/Utilities/Python/fdsplotlib.py b/Utilities/Python/fdsplotlib.py index 769ef2674e5..580e46915f7 100644 --- a/Utilities/Python/fdsplotlib.py +++ b/Utilities/Python/fdsplotlib.py @@ -577,7 +577,7 @@ def plot_to_fig(x_data,y_data,**kwargs): linewidth = kwargs.get('linewidth',1) markeredgewidth = kwargs.get('markeredgewidth',1) markersize = kwargs.get('markersize',5) - + # adjust ticks if required xnumticks = kwargs.get('xnumticks',None) ynumticks = kwargs.get('ynumticks',None) @@ -588,7 +588,7 @@ def plot_to_fig(x_data,y_data,**kwargs): markevery = kwargs.get('data_markevery',default_markevery) legend_location = kwargs.get('legend_location',default_legend_location) legend_framealpha = kwargs.get('legend_framealpha',default_legend_framealpha) - + # set dashes to default, or user requested # This set is the matplotlib default if linestyle == '': dashes = (None, None); linewidth = 0; @@ -596,14 +596,14 @@ def plot_to_fig(x_data,y_data,**kwargs): if linestyle == '--': dashes = kwargs.get('line_dashes',(6, 6)) if linestyle == '-.': dashes = kwargs.get('line_dashes',(6, 3, 1, 3)) if linestyle == ':': dashes = kwargs.get('line_dashes',(1, 3)) - + # This set is what we were using in Matlab # if linestyle == '': dashes = (None, None); linewidth = 0; # if linestyle == '-': dashes = (None, None) # if linestyle == '--': dashes = kwargs.get('line_dashes',(10, 6.2)) # if linestyle == '-.': dashes = kwargs.get('line_dashes',(12, 7.4, 3, 7.4)) # if linestyle == ':': dashes = kwargs.get('line_dashes',(1, 3)) - + data_label = kwargs.get('data_label',None) @@ -695,6 +695,7 @@ def plot_to_fig(x_data,y_data,**kwargs): xerr = kwargs.get('x_error', None) yerr = kwargs.get('y_error', None) + err_linewidth=kwargs.get('error_linewidth', 1) if xerr is not None or yerr is not None: ax.errorbar( x_data, y_data, @@ -704,8 +705,7 @@ def plot_to_fig(x_data,y_data,**kwargs): markeredgewidth=markeredgewidth, # marker edge width markerfacecolor=markerfacecolor, # make marker hollow markeredgecolor=color, # outline color - linestyle=linestyle, - linewidth=linewidth, + elinewidth=err_linewidth, capsize=kwargs.get('error_capsize', 5), # size of caps at ends capthick=linewidth, ) @@ -790,9 +790,9 @@ def plot_to_fig(x_data,y_data,**kwargs): ax.set_xlim(xmin,xmax) ax.set_ylim(ymin,ymax) - + # set number of ticks if requested by the user - if ynumticks != None: + if xnumticks != None: if plot_type in ('loglog', 'semilogx'): ax.set_xticks(np.logspace(xmin, xmax, xnumticks)) else: @@ -802,11 +802,11 @@ def plot_to_fig(x_data,y_data,**kwargs): ax.set_yticks(np.logspace(ymin, ymax, ynumticks)) else: ax.set_yticks(np.linspace(ymin, ymax, ynumticks)) - + # set ticks if requested by the user if xticks is not None: ax.set_xticks(xticks) if yticks is not None: ax.set_yticks(yticks) - + if plot_type in ('loglog', 'semilogy'): apply_global_exponent(ax, axis='y', fontsize=axeslabel_fontsize) @@ -817,7 +817,7 @@ def plot_to_fig(x_data,y_data,**kwargs): add_version_string(ax=ax, version_str=kwargs.get('revision_label'), plot_type=plot_type, font_size=version_fontsize) # fig.tight_layout() # this should not be needed if figure_size and plot_size are both specified - + set_ticks_like_matlab(fig) return fig @@ -1406,6 +1406,9 @@ def scatplot(saved_data, drange, **kwargs): Output_Histograms = [] for _, row in Q.iterrows(): + plt.close('all') + plt.figure().clear() + Scatter_Plot_Title = row["Scatter_Plot_Title"] Plot_Filename = row["Plot_Filename"] Plot_Min = float(row["Plot_Min"]) @@ -1470,7 +1473,9 @@ def scatplot(saved_data, drange, **kwargs): plot_type=Plot_Type, x_min=Plot_Min, x_max=Plot_Max, y_min=Plot_Min, y_max=Plot_Max, x_label=row["Ind_Title"], - y_label=row["Dep_Title"]) + y_label=row["Dep_Title"], + legend_location='outside', + legend_expand=row["Paper_Width_Factor"]) fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=[Plot_Min, Plot_Max], line_style="k-", figure_handle=fig) if Sigma_E > 0: @@ -1480,13 +1485,48 @@ def scatplot(saved_data, drange, **kwargs): fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta * (1 + 2 * Sigma_M), line_style="r--", figure_handle=fig) fdsplotlib.plot_to_fig(x_data=[Plot_Min, Plot_Max], y_data=np.array([Plot_Min, Plot_Max]) * delta / (1 + 2 * Sigma_M), line_style="r--", figure_handle=fig) - # plot data last so it shows on top of stat lines - fdsplotlib.plot_to_fig(x_data=Measured_Values, y_data=Predicted_Values, marker_style="ko", figure_handle=fig) + # --- Plot each dataset with its own marker and color --- + seen_labels = set() + + for idx in match_idx: + style = str(Save_Group_Style[idx]).strip() if Save_Group_Style[idx] else "ko" + fill = str(Save_Fill_Color[idx]).strip() if Save_Fill_Color[idx] else "none" + label = str(Save_Group_Key_Label[idx]).strip() if Save_Group_Key_Label[idx] else "" + + # Flatten valid points for this dataset + mvals = np.array(Save_Measured_Metric[idx], dtype=float).flatten() + pvals = np.array(Save_Predicted_Metric[idx], dtype=float).flatten() + + mask = ( + (mvals >= Plot_Min) & (mvals <= Plot_Max) & + (pvals >= Plot_Min) & (pvals <= Plot_Max) & + (mvals > 0) & (pvals > 0) + ) + mvals = mvals[mask] + pvals = pvals[mask] + + if len(mvals) == 0: + continue + + # Only assign a legend label once per experiment + data_label = label if label and label not in seen_labels else None + if label: + seen_labels.add(label) + + fdsplotlib.plot_to_fig( + x_data=mvals, + y_data=pvals, + marker_style=style, + marker_facecolor=fill, + figure_handle=fig, + data_label=data_label, + ) pdf_path = os.path.join(Manuals_Dir, Plot_Filename + ".pdf") os.makedirs(os.path.dirname(pdf_path), exist_ok=True) fig.savefig(pdf_path) plt.close(fig) + plt.figure().clear() # --- Collect statistics for CSV/TeX --- group_labels = [] @@ -1727,17 +1767,20 @@ def statistics_histogram(Measured_Values, Predicted_Values, outpath = os.path.join(Manuals_Dir, f"{Plot_Filename}_Histogram.pdf") os.makedirs(os.path.dirname(outpath), exist_ok=True) + plt.tight_layout() fig.savefig(outpath) plt.close(fig) + plt.figure().clear() print(f"[statistics_histogram] {Scatter_Plot_Title}: p = {pval:.2f}") return f"{os.path.basename(Plot_Filename)}_Histogram", pval + def set_ticks_like_matlab(fig): ax = fig.axes[0] ax.tick_params(axis="both", direction="in", top=True, right=True, width=0.5) for axis in ['top','bottom','left','right']: ax.spines[axis].set_linewidth(0.5) - + diff --git a/Utilities/Python/scripts/Backward_Facing_Step.py b/Utilities/Python/scripts/Backward_Facing_Step.py new file mode 100644 index 00000000000..80ab2fafa36 --- /dev/null +++ b/Utilities/Python/scripts/Backward_Facing_Step.py @@ -0,0 +1,542 @@ +# Toms +# 8-8-14 +# backward_facing_step.m +# Converted by Floyd +# +# 10/20/2025 + +import numpy as np +import pandas as pd +import os +import math +import matplotlib.pyplot as plt +import fdsplotlib + + +def tight_subplot(Nh, Nw, gap=(0.01, 0.01), marg_h=(0.13, 0.08), marg_w=(0.11, 0.019)): + ''' + Creates a grid of subplot axes with adjustable gaps and margins. + Returns the Figure object and a flat list of Axes handles. + ''' + + if np.isscalar(gap): + gap = (gap, gap) + if np.isscalar(marg_w): + marg_w = (marg_w, marg_w) + if np.isscalar(marg_h): + marg_h = (marg_h, marg_h) + + axh = (1.0 - sum(marg_h) - (Nh - 1) * gap[0]) / Nh + axw = (1.0 - sum(marg_w) - (Nw - 1) * gap[1]) / Nw + + py = 1.0 - marg_h[1] - axh + + ha = [] + fig = plt.figure(figsize=(12, 8)) + + for ih in range(Nh): + px = marg_w[0] + for ix in range(Nw): + # Position: [left, bottom, width, height] + pos = [px, py, axw, axh] + ax = fig.add_axes(pos) + ha.append(ax) + px = px + axw + gap[1] + py = py - axh - gap[0] + + return fig, ha + +expdir = '../../../exp/Backward_Facing_Step/' +datdir = '../../../out/Backward_Facing_Step/' +pltdir = '../../Manuals/FDS_Validation_Guide/SCRIPT_FIGURES/Backward_Facing_Step/' +gitname= 'backward_facing_step_20_git.txt' +git_file = datdir+gitname +version_string = fdsplotlib.get_version_string(git_file) + + +rkappa = 1/.41 +B = 5.2 +rho = 1.1992662 # density, kg/m^3 +U_0 = 7.72 # reference velocity, m/s +h = 0.0098 # step height, m +visc = 1.7801E-5 +visc_nu = visc/rho +x_loc = ['-3','4','6','10'] + +nx = [5, 10, 20] +lnx = len(nx) +fds_marker = ['rs-','go-','bd-'] +subcol = ['r','g','b'] +subsym = ['s','o','d'] +fds_key = ['$h/\\delta x$=5','$h/\\delta x$=10','$h/\\delta x$=20'] + +# --- File Existence Check --- +exp_file = os.path.join(expdir, 'backward_facing_step_data.csv') +if not os.path.exists(exp_file): + print(f'Error: File {exp_file} does not exist. Skipping case.') + quit() + +data_files = [] +for n in nx: + dat_file = os.path.join(datdir, f'backward_facing_step_{n}_line.csv') + if not os.path.exists(dat_file): + print(f'Error: File {dat_file} does not exist. Skipping case.') + quit() + data_files.append(dat_file) + +# --- Read Experimental and FDS Data --- + +# Experimental Data (D): Assumed header is the first row (header=0) +D = pd.read_csv(exp_file) + +# FDS Data (M): Assumed header is the second row (header=1) +M = [] +for df_path in data_files: + M.append(pd.read_csv(df_path, skiprows=1)) + + +# --- 1. Streamwise Cf data along bottom of channel --- + +# Exp Data (Cf-x/h) +xoh_cf = D['Cf-x/h'].values +Cf_exp = D['Cf'].values + +error = 0.0005 * np.ones_like(Cf_exp) + +fig = fdsplotlib.plot_to_fig(x_data=xoh_cf, y_data=Cf_exp, y_error=error,marker_style='ko', + revision_label=version_string,x_min=0,x_max=20,y_min=-0.004,y_max=0.004, + data_label='J&D', + x_label=r'$x/h$', + y_label=r'$C_f$') + +# FDS Data +h_leg = [] +for q in range(lnx): + df = M[q] + + # Column names assumed from MATLAB script: 'u_tau-x' and 'u_tau' + x = df['u_tau-x'].values + I = x > 0 + x = x[I] + + u_tau = df['u_tau'].values[I] + # tau_w = rho * u_tau**2 # Not used for Cf calculation here + + # Column name assumed from MATLAB script: 'u_wall' + u_wall = df['u_wall'].values[I] + + dx_q = h / nx[q] + Cf_fds = visc * u_wall / (0.5 * dx_q) / (0.5 * rho * U_0**2) + + fdsplotlib.plot_to_fig(x_data=x/h, y_data= Cf_fds, marker_style=fds_marker[q],data_markevery=int(len(Cf_fds)/(14+q)), + figure_handle=fig, + data_label=fds_key[q]) + +plotname = pltdir + 'backward_facing_step_Cf.pdf' +plt.savefig(plotname, format='pdf') +plt.close() + +# --- 2. Pressure coefficient along bottom of channel (Cp) --- + +# Exp Data (Cp-x/h) +xoh_cp = D['Cp-x/h'].values +Cp_exp = D['Cp'].values + +# Filter data where x/h <= 20 +I = xoh_cp <= 20 +xoh_cp = xoh_cp[I] +Cp_exp = Cp_exp[I] + +fig = fdsplotlib.plot_to_fig(x_data=xoh_cp, y_data=Cp_exp, marker_style='ko', + revision_label=version_string,x_min=0,x_max=20,y_min=-0.1,y_max=0.25, + data_label='J&D', + x_label=r'$x/h$', + y_label=r'$C_p$') + +# FDS Data +h_leg2 = [] +for q in range(lnx): + df = M[q] + + x = df['cp-x'].values + I = x > 0 + x = x[I] + + Cp_fds_raw = df['cp'].values[I] + + # Cp_fds = Cp / U_0^2 + Cp_fds = Cp_fds_raw / (U_0**2) + + # Normalize Cp data to match experimental end point + # Cp_fds = Cp_fds + Cp(end) - Cp_fds(end); + Cp_fds = Cp_fds + Cp_exp[-1] - Cp_fds[-1] + + + fdsplotlib.plot_to_fig(x_data=x/h, y_data= Cp_fds, marker_style=fds_marker[q],data_markevery=int(len(Cp_fds)/(14+q)), + figure_handle=fig, + data_label=fds_key[q]) + +plotname = pltdir + 'backward_facing_step_Cp.pdf' +plt.savefig(plotname, format='pdf') +plt.close() + +# --- 3. Wall-normal streamwise velocity profiles (U) --- + +fig, sp1 = tight_subplot(1, 4, gap=(0.01, 0.01), marg_h=(0.142, 0.055), marg_w=(0.108, 0.01)) +fig.suptitle(version_string,x=0.8,y=0.98,fontsize=14) + +for i in range(4): # Loop over x_loc: -3, 4, 6, 10 + ax = sp1[i] + x_l = x_loc[i] + + # Exp Data + z_col = f'z {x_l}' + u_col = f'U {x_l}' + + z_exp = D[z_col].values + u_data = D[u_col].values + + I = z_exp > 0 + z_exp = z_exp[I] * 0.001 # Convert mm to m + u_data = u_data[I] + + if i == 0: # x/h = -3 + z_exp += h # Add step height to z coordinate + + ax.plot(u_data / U_0, z_exp / h, 'ko', markersize=10) + + # FDS Data + for q in range(lnx): + df = M[q] + + z_fds_col = f'U-VEL {x_l}-z' + u_fds_col = f'U-VEL {x_l}' + + z_fds = df[z_fds_col].values + u_vel = df[u_fds_col].values + + I = z_fds > 0 + z_fds = z_fds[I] + u_vel = u_vel[I] + + # Plot line and markers + ax.plot( + u_vel/ U_0, + z_fds/ h, + marker=subsym[q], + markevery=int(len(z_fds) /(14+q)), + linestyle='-', + color=subcol[q], + linewidth=1, + markersize=10 + ) + + ax.set_xlim(-0.2, 1.0) + ax.set_ylim(0, 3.5) + ax.tick_params(axis='both', which='major', labelsize=14) + ax.text(0.01, 3.2, f'$x/h$={x_l}', transform=ax.transData, fontsize=14) + + if i == 0: + ax.set_ylabel(r'$z/h$', fontsize=16) + ax.set_xlabel(r'$\langle u \rangle/U_0$', fontsize=16) + # Add legend only to the first plot in this set + legend_labels = ['J&D'] + for k in fds_key: + legend_labels.append(k) + ax.legend(legend_labels, loc='best', frameon=False, fontsize=10) + else: + ax.set_yticklabels([]) + ax.set_xticklabels([]) + +fig.subplots_adjust(left=0.108, right=1.0 - 0.01, bottom=0.055, top=1.0 - 0.142) +fig.savefig(os.path.join(pltdir, 'backward_facing_step_U.pdf')) +plt.close(fig) + +# --- 4. Wall-normal vertical velocity profiles (W) --- + +fig, sp1 = tight_subplot(1, 4, gap=(0.01, 0.01), marg_h=(0.142, 0.055), marg_w=(0.108, 0.01)) +fig.suptitle(version_string,x=0.8,y=0.98,fontsize=14) + +for i in range(4): + ax = sp1[i] + x_l = x_loc[i] + + # Exp Data + z_col = f'z {x_l}' + w_col = f'W {x_l}' + + z_exp = D[z_col].values + w_data = D[w_col].values + + I = z_exp > 0 + z_exp = z_exp[I] * 0.001 + w_data = w_data[I] + + if i == 0: + z_exp += h + + ax.plot(w_data / U_0, z_exp / h, 'ko', markersize=10) + + # FDS Data + for q in range(lnx): + df = M[q] + + z_fds_col = f'W-VEL {x_l}-z' + w_fds_col = f'W-VEL {x_l}' + + z_fds = df[z_fds_col].values + w_vel = df[w_fds_col].values + + I = z_fds > 0 + z_fds = z_fds[I] + w_vel = w_vel[I] + + ax.plot( + w_vel/ U_0, + z_fds/ h, + marker=subsym[q], + markevery=int(len(z_fds) /(14+q)), + linestyle='-', + color=subcol[q], + linewidth=1, + markersize=10 + ) + + ax.set_xlim(-0.1, 0.1) + ax.set_ylim(0, 3.5) + ax.tick_params(axis='both', which='major', labelsize=14) + ax.text(0.01, 3.2, f'$x/h$={x_l}', transform=ax.transData, fontsize=14) + + if i == 0: + ax.set_ylabel(r'$z/h$', fontsize=16) + ax.set_xlabel(r'$\langle w \rangle/U_0$', fontsize=16) + # Add legend only to the first plot in this set + legend_labels = ['J&D'] + for k in fds_key: + legend_labels.append(k) + ax.legend(legend_labels, loc='best', frameon=False, fontsize=10) + else: + ax.set_yticklabels([]) + ax.set_xticklabels([]) + +fig.subplots_adjust(left=0.108, right=1.0 - 0.01, bottom=0.055, top=1.0 - 0.142) +fig.savefig(os.path.join(pltdir, 'backward_facing_step_W.pdf')) +plt.close(fig) + +# --- 5. Wall-normal uu profiles ($\langle uu \rangle$) --- + +fig, sp1 = tight_subplot(1, 4, gap=(0.01, 0.01), marg_h=(0.142, 0.055), marg_w=(0.108, 0.01)) +fig.suptitle(version_string,x=0.8,y=0.98,fontsize=14) + +for i in range(4): + ax = sp1[i] + x_l = x_loc[i] + + # Exp Data + z_col = f'z {x_l}' + uu_col = f'uu {x_l}' + + z_exp = D[z_col].values + uu_data = D[uu_col].values + + I = z_exp > 0 + z_exp = z_exp[I] * 0.001 + uu_data = uu_data[I] + + if i == 0: + z_exp += h + + # No error bars plotted in the original (herrorbar was commented out) + ax.plot(uu_data / U_0**2, z_exp / h, 'ko', markersize=10) + + # FDS Data + for q in range(lnx): + df = M[q] + + z_fds_col = f'uu {x_l}-z' + uu_fds_col = f'uu {x_l}' + + z_fds = df[z_fds_col].values + uu_fds = df[uu_fds_col].values + + I = z_fds > 0 + z_fds = z_fds[I] + uu_fds = uu_fds[I] + + ax.plot( + uu_fds/ U_0**2, + z_fds/ h, + marker=subsym[q], + markevery=int(len(z_fds) /(14+q)), + linestyle='-', + color=subcol[q], + linewidth=1, + markersize=10 + ) + + ax.set_xlim(0, 0.04) + ax.set_ylim(0, 3.5) + ax.tick_params(axis='both', which='major', labelsize=14) + ax.text(0.01, 3.2, f'$x/h={x_l}$', transform=ax.transData, fontsize=14) + + if i == 0: + ax.set_ylabel(r'$z/h$', fontsize=16) + ax.set_xlabel(r'$\langle uu \rangle/U_0^2$', fontsize=16) + # Add legend only to the first plot in this set + legend_labels = ['J&D'] + for k in fds_key: + legend_labels.append(k) + ax.legend(legend_labels, loc='best', frameon=False, fontsize=10) + else: + ax.set_yticklabels([]) + ax.set_xticklabels([]) + +fig.subplots_adjust(left=0.108, right=1.0 - 0.01, bottom=0.055, top=1.0 - 0.142) +fig.savefig(os.path.join(pltdir, 'backward_facing_step_uu.pdf')) +plt.close(fig) + +# --- 6. Wall-normal ww profiles ($\langle ww \rangle$) --- + +fig, sp1 = tight_subplot(1, 4, gap=(0.01, 0.01), marg_h=(0.142, 0.055), marg_w=(0.108, 0.01)) +fig.suptitle(version_string,x=0.8,y=0.98,fontsize=14) + +for i in range(4): + ax = sp1[i] + x_l = x_loc[i] + + # Exp Data + z_col = f'z {x_l}' + ww_col = f'ww {x_l}' + + z_exp = D[z_col].values + ww_data = D[ww_col].values + + I = z_exp > 0 + z_exp = z_exp[I] * 0.001 + ww_data = ww_data[I] + + if i == 0: + z_exp += h + + ax.plot(ww_data / U_0**2, z_exp / h, 'ko', markersize=10) + + # FDS Data + for q in range(lnx): + df = M[q] + + z_fds_col = f'ww {x_l}-z' + ww_fds_col = f'ww {x_l}' + + z_fds = df[z_fds_col].values + ww_fds = df[ww_fds_col].values + + I = z_fds > 0 + z_fds = z_fds[I] + ww_fds = ww_fds[I] + + ax.plot( + ww_fds/ U_0**2, + z_fds/ h, + marker=subsym[q], + markevery=int(len(z_fds) /(14+q)), + linestyle='-', + color=subcol[q], + linewidth=1, + markersize=10 + ) + + ax.set_xlim(0, 0.04) + ax.set_ylim(0, 3.5) + ax.tick_params(axis='both', which='major', labelsize=14) + ax.text(0.01, 3.2, f'$x/h$={x_l}', transform=ax.transData, fontsize=14) + + if i == 0: + ax.set_ylabel(r'$z/h$', fontsize=16) + ax.set_xlabel(r'$\langle ww \rangle/U_0^2$', fontsize=16) + # Add legend only to the first plot in this set + legend_labels = ['J&D'] + for k in fds_key: + legend_labels.append(k) + ax.legend(legend_labels, loc='best', frameon=False, fontsize=10) + else: + ax.set_yticklabels([]) + ax.set_xticklabels([]) + +fig.subplots_adjust(left=0.108, right=1.0 - 0.01, bottom=0.055, top=1.0 - 0.142) +fig.savefig(os.path.join(pltdir, 'backward_facing_step_ww.pdf')) +plt.close(fig) + +fig, sp1 = tight_subplot(1, 4, gap=(0.01, 0.01), marg_h=(0.142, 0.055), marg_w=(0.108, 0.01)) +fig.suptitle(version_string,x=0.8,y=0.98,fontsize=14) + +# --- 7. Wall-normal uw profiles ($\langle uw \rangle$) --- + +for i in range(4): + ax = sp1[i] + x_l = x_loc[i] + + # Exp Data + z_col = f'z {x_l}' + uw_col = f'uw {x_l}' + + z_exp = D[z_col].values + uw_data = D[uw_col].values + + I = z_exp > 0 + z_exp = z_exp[I] * 0.001 + uw_data = uw_data[I] + + if i == 0: + z_exp += h + + ax.plot(uw_data / U_0**2, z_exp / h, 'ko', markersize=10) + + # FDS Data + h_leg = [] + for q in range(lnx): + df = M[q] + + z_fds_col = f'uw {x_l}-z' + uw_fds_col = f'uw {x_l}' + + z_fds = df[z_fds_col].values + uw_fds_raw = df[uw_fds_col].values + + I = z_fds > 0 + z_fds = z_fds[I] + uw_fds_raw = uw_fds_raw[I] + + # Original MATLAB applies -1 multiplication: uw_fds = -1*M{q}.data(I,j); + uw_fds = -1 * uw_fds_raw + ax.plot( + uw_fds/ U_0**2, + z_fds/ h, + marker=subsym[q], + markevery=int(len(z_fds) /(14+q)), + linestyle='-', + color=subcol[q], + linewidth=1, + markersize=10 + ) + + ax.set_xlim(0, 0.04) # Note: Original MATLAB had axis([0 .04 0 3.5]) + ax.set_ylim(0, 3.5) + ax.tick_params(axis='both', which='major', labelsize=14) + ax.text(0.01, 3.2, f'$x/h$={x_l}', transform=ax.transData, fontsize=14) + + if i == 0: + ax.set_ylabel(r'$z/h$', fontsize=16) + ax.set_xlabel(r'$\langle uw \rangle/U_0^2$', fontsize=16) + # Add legend only to the first plot in this set + legend_labels = ['J&D'] + for k in fds_key: + legend_labels.append(k) + ax.legend(legend_labels, loc='best', frameon=False, fontsize=10) + else: + ax.set_yticklabels([]) + ax.set_xticklabels([]) + +fig.subplots_adjust(left=0.108, right=1.0 - 0.01, bottom=0.055, top=1.0 - 0.142) +fig.savefig(os.path.join(pltdir, 'backward_facing_step_uw.pdf')) +plt.close(fig) + diff --git a/Utilities/Python/scripts/atmospheric_boundary_layer.py b/Utilities/Python/scripts/atmospheric_boundary_layer.py index 177bafb32c8..eb7977e220a 100644 --- a/Utilities/Python/scripts/atmospheric_boundary_layer.py +++ b/Utilities/Python/scripts/atmospheric_boundary_layer.py @@ -11,7 +11,7 @@ plot_style = fdsplotlib.get_plot_style('fds') outdir = '../../Verification/Atmospheric_Effects/' -pltdir = '../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/' +pltdir = '../../Manuals/FDS_User_Guide/SCRIPT_FIGURES/' git_file = outdir + 'atmospheric_boundary_layer_1_git.txt' version_string = fdsplotlib.get_version_string(git_file) diff --git a/Utilities/Python/scripts/catchpole_spread_rates.py b/Utilities/Python/scripts/catchpole_spread_rates.py index dc65c5e845b..4b090ba76b5 100644 --- a/Utilities/Python/scripts/catchpole_spread_rates.py +++ b/Utilities/Python/scripts/catchpole_spread_rates.py @@ -1,3 +1,6 @@ +# This script creates additional figures related to the USFS_Catchpole cases. +# The objective is to better identify trends/issues within the parameter space. + import numpy as np import matplotlib.pyplot as plt import os, sys @@ -22,24 +25,21 @@ # Experiment parameters tests = pd.read_csv(os.path.join(validation_path,"Test_Matrix.csv")) +# Calculate spread rates for summary plotting for ti,test in tests.iterrows(): R = tests['R'].iloc[ti] chid = tests['Test'].iloc[ti] fds_file = os.path.join(base_path, f"{chid}_devc.csv") git_file = os.path.join(base_path, f"{chid}_git.txt") - fig_file = os.path.join(fig_path, f"{chid}.pdf") - - if os.path.exists(fds_file) is False: - print(f'Error: File {fds_file} does not exist. Skipping case.') - continue + version_string = fdsplotlib.get_version_string(git_file) fds_data = pd.read_csv(fds_file,header=1) # fit spread rate try: - # fit slope filtering to positions greater than 2 m from ignition - R_FDS,intercept = np.polyfit(fds_data[fds_data['x']>2]['Time'], - fds_data[fds_data['x']>2]['x'], 1) # 1 indicates linear fit (degree 1) + # fit slope filtering to positions between 2 m and 7 m from ignition + xi = (fds_data['x']>2) & (fds_data['x']<7) + R_FDS,intercept = np.polyfit(fds_data[xi]['Time'], fds_data[xi]['x'], 1) # not enough data to fit except: R_FDS=0. @@ -47,30 +47,14 @@ if R_FDS<0: R_FDS=0 - x_exp = np.array([0., 8.]) - t_exp = np.array([0., 8./R]) - version_string = fdsplotlib.get_version_string(git_file) - fig = fdsplotlib.plot_to_fig(x_data=t_exp, y_data=x_exp, - revision_label=version_string, - x_min=0., x_max=8./R, y_min=0., y_max=8., - x_label='Time (s)', y_label='Distance (m)', - data_label='EXP', - marker_style='k-') - - fdsplotlib.plot_to_fig(x_data=fds_data['Time'].values, y_data=fds_data['x'].values, figure_handle=fig, data_label='FDS', marker_style='k--') - - plt.savefig(fig_file) - plt.close() - - # write table for dataplot - tests.loc[tests['Test'].index[ti],'R_FDS'] = R_FDS - test = tests.iloc[ti] - test = test.drop('Test') - out_file = os.path.join(base_path,f"{chid}_FDS.csv") - pd.DataFrame([test]).to_csv(out_file,index=False) - # add fds data to full table for summary plotting tests.loc[ti,'R_FDS'] = R_FDS + # No spread if zero from above or the fire does not reach at least 7 m + if (R_FDS < 1e-5) or (fds_data['x'].max() < 7): + tests.loc[ti,'category'] = 'no spread' + else: + tests.loc[ti,'category'] = 'spread' + # Create summary plots @@ -86,7 +70,7 @@ fig_file = os.path.join(fig_path, f"Catchpole_R_v_{dvar}.pdf") # show +/- 20% relative error - [xmin,xmax] = [tests[dvar].min(),tests[dvar].max()] + [xmin,xmax] = [0.8*tests[dvar].min(),1.1*tests[dvar].max()] fig = fdsplotlib.plot_to_fig(x_data=[xmin, xmax], y_data=[0.8,0.8], plot_type='semilogy', marker_style='k--', x_min=xmin, x_max=xmax, y_min=3e-3, y_max=1.1e1, @@ -98,7 +82,6 @@ for i in range(0, len(fuel_labels)): fuel = fuel_labels[i] filtered_data = tests[tests['Test'].str.startswith(fuel)] - color = colors[i] if fuel=='EX': filtered_data = tests[ (tests['Test'].str.startswith(fuel))&(~tests['Test'].str.startswith('EXSC'))] @@ -111,29 +94,25 @@ # plot no-spread conditions fig_file = os.path.join(fig_path, "Catchpole_no_spread.pdf") -fig, ax = plt.subplots(figsize=(plot_style["Paper_Width"], plot_style["Paper_Height"])) - -# dummy column for labeling -tests['category']='go' -tests.loc[tests['R_FDS']<1e-5,'category'] = 'no-go' - -# normalize by max and min +# Dummy call to establish figure +fig = fdsplotlib.plot_to_fig(x_data=[-1,-1], y_data=[-1,-1], + x_label='Parameter',y_label='Normalized value', + x_min=0,x_max=3,y_min=0,y_max=1,ynumticks=2, + revision_label=version_string) +ax = plt.gca() + +# Normalize by max and min tests_normalized = tests tests_normalized[list(dep_variables.keys())] = tests[list(dep_variables.keys())].apply( lambda x: (x - x.min()) / (x.max() - x.min())) -# move M toward the middle of x-axis for more clarity -tests_normalized=tests_normalized[['category','s','beta','M','U']] - pd.plotting.parallel_coordinates(tests_normalized, 'category', cols=['s','beta','M','U'], color=[(1.,0.,0.,1), (0.,0.,0.,.2)], ax=ax, ls='-') -ax.set_ylim([0, 1]) -ax.set_yticks([0, 1],['min','max']) -plt.legend(loc="upper left", fontsize=plot_style["Key_Font_Size"], framealpha=1,frameon=True) +ax.legend(loc="upper center", framealpha=1,frameon=True) plt.savefig(fig_file) plt.close() diff --git a/Utilities/Python/scripts/compression_wave.py b/Utilities/Python/scripts/compression_wave.py index 304d464e779..48c32ebed22 100644 --- a/Utilities/Python/scripts/compression_wave.py +++ b/Utilities/Python/scripts/compression_wave.py @@ -17,25 +17,40 @@ def compression_wave_soln(rho0, x, y, a, c, t): + x = np.asarray(x) + y = np.asarray(y) + rho0 = np.asarray(rho0) - b = np.sqrt(-1 + a**2) - d = np.sqrt(-1 + c**2) + b = np.sqrt(a**2 - 1.0) + d = np.sqrt(c**2 - 1.0) - x0 = 2 * np.arctan(b/a * np.tan(np.arctan((1 + a*np.tan(x/2))/b) - b*t/2) - 1/a) - y0 = 2 * np.arctan(d/c * np.tan(np.arctan((1 + c*np.tan(y/2))/d) - d*t/2) - 1/c) + x0 = 2*np.arctan( (b/a)*np.tan(np.arctan((1 + a*np.tan(x/2))/b) - (b*t)/2) - 1/a ) + y0 = 2*np.arctan( (d/c)*np.tan(np.arctan((1 + c*np.tan(y/2))/d) - (d*t)/2) - 1/c ) - Ix0 = np.log(-a**2 - np.cos(2*np.arctan((1 + a*np.tan(x0/2))/b)) + b*np.sin(2*np.arctan((1 + a*np.tan(x0/2))/b))) - Iy0 = np.log(-c**2 - np.cos(2*np.arctan((1 + c*np.tan(y0/2))/d)) + d*np.sin(2*np.arctan((1 + c*np.tan(y0/2))/d))) + # small positive to avoid exact zeros inside log + tiny = np.finfo(float).tiny - Ix = np.log(-a**2 - np.cos(b*t + 2*np.arctan((1 + a*np.tan(x0/2))/b)) + b*np.sin(b*t + 2*np.arctan((1 + a*np.tan(x0/2))/b))) - Iy = np.log(-c**2 - np.cos(d*t + 2*np.arctan((1 + c*np.tan(y0/2))/d)) + d*np.sin(d*t + 2*np.arctan((1 + c*np.tan(y0/2))/d))) + def clog(real_expr): + # promote to complex; add tiny where exactly zero to avoid log(0) + z = real_expr.astype(np.complex128) + z = np.where(real_expr == 0, tiny + 0j, z) + return np.log(z) - q0 = np.log(rho0) - q = q0 + Ix - Ix0 + Iy - Iy0 + A0 = 2*np.arctan((1 + a*np.tan(x0/2))/b) + C0 = 2*np.arctan((1 + c*np.tan(y0/2))/d) + + Ix0 = clog(-a**2 - np.cos(A0) + b*np.sin(A0)) + Iy0 = clog(-c**2 - np.cos(C0) + d*np.sin(C0)) - rho = np.exp(q) + Ix = clog(-a**2 - np.cos(b*t + A0) + b*np.sin(b*t + A0)) + Iy = clog(-c**2 - np.cos(d*t + C0) + d*np.sin(d*t + C0)) + + q0 = np.log(rho0) + q = q0 + (Ix - Ix0) + (Iy - Iy0) - return rho + rho = np.exp(q) # complex in general + rho = np.real_if_close(rho, tol=1e6) # drop ~0 imag parts (e.g., ≤1e-10) + return rho.real # central differencing, FL=0 diff --git a/Utilities/Python/scripts/geom_channel_test.py b/Utilities/Python/scripts/geom_channel_test.py new file mode 100644 index 00000000000..17f7fae9ab9 --- /dev/null +++ b/Utilities/Python/scripts/geom_channel_test.py @@ -0,0 +1,59 @@ +#!/usr/bin/env python +""" +geom_channel_test.py +Temporary home for checking AREA integration for geom_channel.fds +Original MATLAB script by McDermott (8-11-2021) +""" +import numpy as np +import pandas as pd +import os + +def check_geom_channel(): + """Check AREA integration for geom_channel test case""" + + fds_dir = os.path.normpath(os.path.join(os.path.dirname(__file__),'..','..','..')) + ddir = os.path.join(fds_dir, 'Verification','Complex_Geometry','') + devc_file = os.path.join(ddir, 'geom_channel_devc.csv') + try: + M = pd.read_csv(devc_file, skiprows=1) # Skip units row, use column names as header + except FileNotFoundError: + print(f'Error: File {devc_file} does not exist. Skipping case.') + return False + + U0 = M['U0'].iloc[-1] + A0 = M['A0'].iloc[-1] + U1 = M['U1'].iloc[-1] + A1 = M['A1'].iloc[-1] + U2 = M['U2'].iloc[-1] + A2 = M['A2'].iloc[-1] + U3 = M['U3'].iloc[-1] + A3 = M['A3'].iloc[-1] + + error_tolerance = 1e-6 + all_passed = True + A0_error = abs(A0 - 1) + if A0_error > error_tolerance: + print(f'Python Warning: geom_channel.fds A0_error = {A0_error:.6e}') + all_passed = False + A1_error = abs(A1 - 1) + if A1_error > error_tolerance: + print(f'Python Warning: geom_channel.fds A1_error = {A1_error:.6e}') + all_passed = False + A2_error = abs(A2 - 1) + if A2_error > error_tolerance: + print(f'Python Warning: geom_channel.fds A2_error = {A2_error:.6e}') + all_passed = False + A3_error = abs(A3 - 1) + if A3_error > error_tolerance: + print(f'Python Warning: geom_channel.fds A3_error = {A3_error:.6e}') + all_passed = False + return all_passed + +if __name__ == '__main__': + """Main execution block""" + print("Running geom channel test verification...") + success = check_geom_channel() + if success: + print("Geom channel test verification completed successfully!") + else: + print("Geom channel test verification completed with errors.") \ No newline at end of file diff --git a/Utilities/Python/scripts/geom_positive_errors.py b/Utilities/Python/scripts/geom_positive_errors.py new file mode 100644 index 00000000000..86de8e8022b --- /dev/null +++ b/Utilities/Python/scripts/geom_positive_errors.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python +""" +geom_positive_errors.py +This script catches errors produced during FDS setup of bad geometries cases. +The success condition is the existence of the error message. The error message has +its 'ERROR' label replaced by 'SUCCESS' label, thanks to the MISC parameter +POSITIVE_ERROR_TEST=.TRUE. in the bad geometries cases. +The absence of the 'ERROR' label makes current Firebot to consider the run a success. +Original MATLAB script by Emanuele Gissi (6-7-2017) +""" +import os + +def check_positive_errors(): + """Check for expected error messages in bad geometry test case error files""" + fds_dir = os.path.normpath(os.path.join(os.path.dirname(__file__),'..','..','..')) + err_dir = os.path.join(fds_dir, 'Verification', 'Complex_Geometry', '') + cases = [ + { + 'file': 'geom_bad_inconsistent_normals.err', + 'error_string': 'Face normals are probably pointing in the wrong direction' + }, + { + 'file': 'geom_bad_open_surface.err', + 'error_string': 'Open geometry at edge' + }, + { + 'file': 'geom_bad_non_manifold_edge.err', + 'error_string': 'Non manifold geometry in adjacent faces at edge' + }, + { + 'file': 'geom_bad_inverted_normals.err', + 'error_string': 'Face normals are probably pointing in the wrong direction' + }, + { + 'file': 'geom_bad_non_manifold_vert.err', + 'error_string': 'Vertex 6 not connected.' + } + ] + all_passed = True + + # Check existence of error_string in each file + for n, case in enumerate(cases, start=1): + infile = os.path.join(err_dir, case['file']) + errstring = case['error_string'] + + # Check if file exists + if not os.path.exists(infile): + print(f'Error: File {infile} does not exist. Skipping case.') + all_passed = False + continue + + # Read file and check for error string + try: + with open(infile, 'r') as f: + errfile_content = f.read() + + if errstring not in errfile_content: + print(f'Error: File {infile}:') + print(f' Does not contain the following positive error message:') + print(f' <{errstring}>') + all_passed = False + except Exception as e: + print(f'Error reading file {infile}: {e}') + all_passed = False + return all_passed + +if __name__ == '__main__': + """Main execution block""" + print("Running geometry positive errors verification...") + success = check_positive_errors() + if success: + print("Geometry positive errors verification completed successfully!") + else: + print("Geometry positive errors verification completed with errors.") \ No newline at end of file diff --git a/Utilities/Python/scripts/mass_balance.py b/Utilities/Python/scripts/mass_balance.py index ca316a41217..7731c727d07 100644 --- a/Utilities/Python/scripts/mass_balance.py +++ b/Utilities/Python/scripts/mass_balance.py @@ -13,8 +13,9 @@ def plot_mass_balance(chid, title_text): plot_style = fdsplotlib.get_plot_style('fds') - ddir = '../../Verification/Species/' - pltdir = '../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/' + fds_dir = os.path.normpath(os.path.join(os.path.dirname(__file__),'..','..','..')) + ddir = os.path.join(fds_dir, 'Verification','Species','') + pltdir = os.path.join(fds_dir, 'Manuals','FDS_Verification_Guide','SCRIPT_FIGURES','') mass_file = os.path.join(ddir, f'{chid}_mass.csv') @@ -47,16 +48,19 @@ def plot_mass_balance(chid, title_text): # Plot zero reference line first (black solid line, no label) fig = fdsplotlib.plot_to_fig(x_data=t, y_data=np.zeros_like(t), marker_style='k-',revision_label=version_string,legend_location='upper right', - x_label='time (s)', y_label='mass flow (kg/s)',x_min=0, x_max=2000, y_min=-0.005, y_max=0.015) + x_label='time (s)', y_label='mass flow (kg/s)',x_min=0, x_max=2000, y_min=-0.005, y_max=0.015, + xnumticks=6, ynumticks=5) fdsplotlib.plot_to_fig(x_data=t, y_data=mdot_in, figure_handle=fig, marker_style='r-', data_label='Inlet H2O') fdsplotlib.plot_to_fig(x_data=t, y_data=-mdot_out, figure_handle=fig, marker_style='m-', data_label='Outlet H2O') fdsplotlib.plot_to_fig(x_data=t, y_data=bal, figure_handle=fig, marker_style='b-', data_label='dm/dt+out-in') ax = plt.gca() lines = ax.get_lines() - lines[2].set_color('#EDB120') # Orange for outlet (second data line) - ax.locator_params(axis='x', nbins=6) # ~6 ticks on x-axis - ax.locator_params(axis='y', nbins=5) # ~5 ticks on y-axis + lines[2].set_color('orange') # Orange for outlet (second data line) + legend = ax.get_legend() + if legend: + legend.legend_handles[1].set_color('orange') # Second legend entry (Outlet H2O) + fdsplotlib.apply_global_exponent(ax, axis='y') ax.text(100, 18e-3, title_text,fontsize=plot_style['Title_Font_Size'],fontname=plot_style['Font_Name']) diff --git a/Utilities/Python/scripts/mass_balance_gas_volume.py b/Utilities/Python/scripts/mass_balance_gas_volume.py new file mode 100644 index 00000000000..271012c4ee8 --- /dev/null +++ b/Utilities/Python/scripts/mass_balance_gas_volume.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python +""" +mass_balance_gas_volume.py +Validates mass conservation in a gas volume comparing time derivative of mass with net mass flux through boundaries. +Original MATLAB script by McDermott (2 Dec 2021) +""" +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import fdsplotlib +import os + +def plot_mass_balance(chid, title_text): + + plot_style = fdsplotlib.get_plot_style('fds') + fds_dir = os.path.normpath(os.path.join(os.path.dirname(__file__),'..','..','..')) + ddir = os.path.join(fds_dir, 'Verification','Species','') + pltdir = os.path.join(fds_dir, 'Manuals','FDS_Verification_Guide','SCRIPT_FIGURES','') + devc_file = os.path.join(ddir, f'{chid}_devc.csv') + + try: + M = pd.read_csv(devc_file, skiprows=1) # Skip units row, use column names as header + except FileNotFoundError: + print(f'Error: File {devc_file} does not exist. Skipping case.') + return 0 + + t = M.iloc[:, 0].values + m = M['M'].values + dmdt = np.zeros(len(t)) + for i in range(1, len(t)): + dmdt[i] = (m[i] - m[i-1]) / (t[i] - t[i-1]) + + mf_x1 = M['TMF_X1'].values + mf_x2 = M['TMF_X2'].values + mf_y1 = M['TMF_Y1'].values + mf_y2 = M['TMF_Y2'].values + mf_z1 = M['TMF_Z1'].values + mf_z2 = M['TMF_Z2'].values + # Net flux: out - in + net_flux = mf_x2 - mf_x1 + mf_y2 - mf_y1 + mf_z2 - mf_z1 + # Residual: dm/dt + (out - in) (should be ~0) + bal = dmdt + net_flux + + # Get version string if git file exists + git_file = os.path.join(ddir, f'{chid}_git.txt') + version_string = fdsplotlib.get_version_string(git_file) if os.path.exists(git_file) else '' + + # Plot zero reference line first (black solid line, no label) + fig = fdsplotlib.plot_to_fig(x_data=t, y_data=np.zeros_like(t), marker_style='k-', + revision_label=version_string, legend_location='upper right',x_label='time (s)', y_label='mass flow (kg/s)', + x_min=0, x_max=2000, y_min=-0.004, y_max=0.004,xnumticks=5, ynumticks=5) + fdsplotlib.plot_to_fig(x_data=t, y_data=dmdt, figure_handle=fig, marker_style='r-', data_label='dm/dt') + fdsplotlib.plot_to_fig(x_data=t, y_data=net_flux, figure_handle=fig, marker_style='m-', data_label='out-in') + fdsplotlib.plot_to_fig(x_data=t, y_data=bal, figure_handle=fig, marker_style='b-', data_label='dm/dt+out-in') + + ax = plt.gca() + lines = ax.get_lines() + lines[2].set_color('orange') # Orange for net flux (second data line) + legend = ax.get_legend() + if legend: + legend.legend_handles[1].set_color('orange') # Second legend entry (out-in) + fdsplotlib.apply_global_exponent(ax, axis='y') + if title_text: + ax.text(100, 18e-3, title_text, fontsize=plot_style['Title_Font_Size'], + fontname=plot_style['Font_Name']) + + # Check balance and report error + mass_error = np.max(np.abs(bal)) + if mass_error > 1e-6: + print(f'Python Warning: mass error = {mass_error:.6e} in {chid}') + + # Save figure + output_file = os.path.join(pltdir, f'{chid}.pdf') + plt.savefig(output_file, format='pdf') + plt.close(fig) + + return mass_error + +if __name__ == '__main__': + """Main execution block""" + print("Running mass balance gas volume verification plot...") + error = plot_mass_balance('mass_balance_gas_volume', '') + print("Mass balance gas volume verification completed successfully!") \ No newline at end of file diff --git a/Utilities/Python/scripts/mass_balance_reac.py b/Utilities/Python/scripts/mass_balance_reac.py new file mode 100644 index 00000000000..93ea79b0b33 --- /dev/null +++ b/Utilities/Python/scripts/mass_balance_reac.py @@ -0,0 +1,87 @@ +#!/usr/bin/env python +""" +mass_balance_reac.py +Validates mass conservation for species with chemical reactions. +Original MATLAB script by McDermott (10 July 2020) +""" +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +import fdsplotlib +import os + +def plot_mass_balance(chid, title_text, mass_id, devc_id, y_min=None, y_max=None, ynumticks=8): + + plot_style = fdsplotlib.get_plot_style('fds') + fds_dir = os.path.normpath(os.path.join(os.path.dirname(__file__),'..','..','..')) + ddir = os.path.join(fds_dir, 'Verification','Species','') + pltdir = os.path.join(fds_dir, 'Manuals','FDS_Verification_Guide','SCRIPT_FIGURES','') + + mass_file = os.path.join(ddir, f'{chid}_mass.csv') + devc_file = os.path.join(ddir, f'{chid}_devc.csv') + + try: + M = pd.read_csv(mass_file, skiprows=1) # Skip units row, use column names as header + F = pd.read_csv(devc_file, skiprows=1) # Skip units row, use column names as header + except FileNotFoundError as e: + print(f'Error: File does not exist. Skipping case. {e}') + return 0 + + t = M.iloc[:, 0].values + m = M[mass_id].values + dmdt = np.zeros(len(t)) + for i in range(1, len(t)): + dmdt[i] = (m[i] - m[i-1]) / (t[i] - t[i-1]) + + mdot_out = (F[f'{devc_id} xmax'].values + + F[f'{devc_id} xmin'].values + + F[f'{devc_id} ymin'].values + + F[f'{devc_id} ymax'].values + + F[f'{devc_id} Burner'].values) + + # Generation/consumption from reaction + gen = F[f'{devc_id} mdot reac'].values + # Balance: -dm/dt + (in - out) + generation (should be ~0) + bal = -dmdt + mdot_out + gen + + # Get version string if git file exists + git_file = os.path.join(ddir, f'{chid}_git.txt') + version_string = fdsplotlib.get_version_string(git_file) if os.path.exists(git_file) else '' + + # Create plot with all data + fig = fdsplotlib.plot_to_fig(x_data=t, y_data=dmdt, marker_style='g-', + revision_label=version_string, legend_location='lower left', x_label='Time (s)', y_label='Mass Flow (kg/s)', + data_label='accumulation',x_min=0, x_max=20, y_min=y_min, y_max=y_max, xnumticks=5, ynumticks=ynumticks) + fdsplotlib.plot_to_fig(x_data=t, y_data=mdot_out, figure_handle=fig, marker_style='b-', data_label='in - out') + fdsplotlib.plot_to_fig(x_data=t, y_data=gen, figure_handle=fig, marker_style='r-', data_label='generation') + fdsplotlib.plot_to_fig(x_data=t, y_data=bal, figure_handle=fig, marker_style='k-', data_label='balance') + ax = plt.gca() + xl = ax.get_xlim() + yl = ax.get_ylim() + x_pos = xl[0] + 0.05 * (xl[1] - xl[0]) + y_pos = yl[0] + 0.9 * (yl[1] - yl[0]) + ax.text(x_pos, y_pos, title_text, fontsize=plot_style['Title_Font_Size'], + fontname=plot_style['Font_Name']) + + # Check balance and report error (normalized by total mass) + mass_error = np.abs(np.mean(bal)) / m[-1] if m[-1] != 0 else np.nan + if mass_error > 1.5e-3: + print(f'Python Warning: mass error = {mass_error:.6e} in {chid} for species {mass_id}') + + # Save figure + output_file = os.path.join(pltdir, f'{chid}_{devc_id}.pdf') + plt.savefig(output_file, format='pdf') + plt.close(fig) + + return mass_error + +if __name__ == '__main__': + """Main execution block""" + print("Running mass balance reac verification plots...") + error1 = plot_mass_balance('mass_balance_reac', 'Propane Mass Balance', 'PROPANE', 'C3H8',y_min=-0.15,y_max=0.1,ynumticks=6) + error2 = plot_mass_balance('mass_balance_reac', 'Oxygen Mass Balance', 'OXYGEN', 'O2',y_min=-2,y_max=1.5,ynumticks=8) + error3 = plot_mass_balance('mass_balance_reac', 'Nitrogen Mass Balance', 'NITROGEN', 'N2',y_min=-6,y_max=4,ynumticks=6) + error4 = plot_mass_balance('mass_balance_reac', 'Carbon Dioxide Mass Balance', 'CARBON DIOXIDE', 'CO2',y_min=-0.4,y_max=0.4,ynumticks=5) + error5 = plot_mass_balance('mass_balance_reac', 'Water Vapor Mass Balance', 'WATER VAPOR', 'H2O',y_min=-0.2,y_max=0.2,ynumticks=5) + error6 = plot_mass_balance('mass_balance_reac', 'Soot Mass Balance', 'SOOT', 'Soot',y_min=-0.002,y_max=0.002,ynumticks=5) + print("Mass balance reac verification completed successfully!") \ No newline at end of file diff --git a/Utilities/Python/scripts/rotcube_cc_mms_error.py b/Utilities/Python/scripts/rotcube_cc_mms_error.py new file mode 100644 index 00000000000..c1f19a83298 --- /dev/null +++ b/Utilities/Python/scripts/rotcube_cc_mms_error.py @@ -0,0 +1,349 @@ +# Vanella rotcube_cc_mms_error.m +# 05-24-2019 +# Computes L1, L2, and Linf errors for the Rotated Cube MMS case. +# +# Converted by Floyd +# 10/17/2025 + +import numpy as np +import pandas as pd +import os +import math +import matplotlib.pyplot as plt +import fdsplotlib + + +mu = 0.01 +D = 0.01 +rho = 1.00 +L = np.pi # Cube side length. + +nwave = 1. # Wave number on analytical solution. +gam = np.pi/2 +Az = 0.1 +meanz = 0.15 +displ = np.pi +dispxy = np.array([-1/2 * displ, -1/2 * displ]) + +datadir = '../../Verification/Complex_Geometry/' +plotdir = '../../Manuals/FDS_Verification_Guide/SCRIPT_FIGURES/' +gitname = 'rotated_cube_45deg_256_stm_git.txt' +git_file = datadir+gitname +version_string = fdsplotlib.get_version_string(git_file) + +def ann_z(x, y, t): + return ( + Az / 3. * np.sin(t) * (1. - np.cos(2. * nwave * (x - gam))) * (1. - np.cos(2. * nwave * (y - gam))) + - Az / 3. * np.sin(t) + + meanz + ) + +def ann_u(x, y, t): + return -np.sin(t) * np.sin(nwave * x)**2 * np.sin(2. * nwave * y) + +def ann_v(x, y, t): + return np.sin(t) * np.sin(2. * nwave * x) * np.sin(nwave * y)**2 + +def ann_p(x, y, t): + return np.sin(t) / 4. * (2. + np.cos(2. * nwave * x)) * (2. + np.cos(2. * nwave * y)) - np.sin(t) + +def ann_H(x, y, t): + u = ann_u(x, y, t) + v = ann_v(x, y, t) + return (u**2 + v**2) / 2 + ann_p(x, y, t) / rho + + +# Function to read the MATLAB-style CSV file header +def read_matlab_csv(filepath): + try: + # Read the first line (metadata) + with open(filepath, 'r') as f: + header_str = f.readline().strip() + + # MATLAB's str2num is flexible; we try to handle space/comma separation + # and ensure we get 6 numbers. + vec = [float(x) for x in header_str.replace(',', ' ').split() if x.strip()] + + if len(vec) < 6: + raise ValueError(f"Could not parse 6 header values from first line: '{header_str}'") + + # Metadata: ntot_u, ntot_v, ntot_c, T, dx, dy + ntot_u, ntot_v, ntot_c = int(vec[0]), int(vec[1]), int(vec[2]) + T, dx, dy = vec[3], vec[4], vec[5] + + M_data = pd.read_csv(filepath, skiprows=1, header=None, usecols=range(7)).values + + return M_data, ntot_u, ntot_v, ntot_c, T, dx, dy + + except FileNotFoundError: + print(f" [ERROR] File not found: {filepath}") + return None, 0, 0, 0, 0, 0, 0 + except Exception as e: + print(f" [ERROR] Processing {filepath}: {e}") + return None, 0, 0, 0, 0, 0, 0 + + +# --- 3. File setup arrays (0-indexed in Python) --- +ifile_s = 1 +ifile_f = 3 +ifile_ang = [0, np.arctan(1/2.), np.arctan(1.)] +ifile_str = ['0deg_', '27deg_', '45deg_'] +jfile_str = ['stm', 'obs'] + +# List of dictionaries to replace MATLAB struct 'file' +files = [] +skip_case = False + +for ifile in range(ifile_s - 1, ifile_f): # indices 0, 1, 2 + jfile_s = 0 + jfile_f = 0 + if ifile == 0: # 0 deg case (stm and obs) + jfile_f = 1 + for jfile in range(jfile_s, jfile_f + 1): # indices 0, 1 (for ifile 0), index 0 (for others) + + # Base file names for resolutions 32, 64, 128, 256, 320 + resolutions = [32, 64, 128, 256, 320] + name_list = [] + for res in resolutions: + name_list.append( + f'rotated_cube_{ifile_str[ifile]}{res}_{jfile_str[jfile]}_mms.csv' + ) + + file_entry = { + 'name': name_list, + 'nameout': f'rotated_cube_{ifile_str[ifile]}{jfile_str[jfile]}_mms_convergence', + 'rotang': ifile_ang[ifile], + 'dx': [], + 'errors': {} # Dictionary to hold all computed error arrays + } + + # Check existence and populate files list + for name in file_entry['name']: + if not os.path.exists(os.path.join(datadir, name)): + skip_case = True + print(f"[WARN] File {os.path.join(datadir, name)} does not exist. Skipping case.") + + files.append(file_entry) + +if skip_case: quit() + +# --- 4. Compute errors from *_mms.csv files --- +for ifile in range(len(files)): + + file_entry = files[ifile] + rotangle = file_entry['rotang'] + + ROTMAT = np.array([ + [np.cos(rotangle), -np.sin(rotangle)], + [np.sin(rotangle), np.cos(rotangle)] + ]) + TROTMAT = ROTMAT.T # Transpose + + # Initialize lists to store errors for all resolutions (n=1 to 5) + error_lists = {k: [] for k in ['e_u_1', 'e_u_2', 'e_u_i', 'e_v_1', 'e_v_2', 'e_v_i', + 'e_z_1', 'e_z_2', 'e_z_i', 'e_H_1', 'e_H_2', 'e_H_i', + 'e_p_1', 'e_p_2', 'e_p_i']} + + for n in range(len(file_entry['name'])): + filepath = os.path.join(datadir, file_entry['name'][n]) + + M_data, ntot_u, ntot_v, ntot_c, T, dx, dy = read_matlab_csv(filepath) + + if M_data is None or M_data.size == 0: + print(f" Skipping processing for {file_entry['name'][n]} due to missing or empty data.") + continue # Skip this resolution + + file_entry['dx'].append(dx) + + # Initialize error vectors + e_u_vec = np.zeros(ntot_u) + e_v_vec = np.zeros(ntot_v) + e_z_vec = np.zeros(ntot_c) + e_H_vec = np.zeros(ntot_c) + e_p_vec = np.zeros(ntot_c) + + p = 0 # Data row counter (0-indexed) + + # --- First U: Velocity component along x-axis --- + # M_data columns: 0: iscf, 1: xglob, 2: yglob, 3: area, 4: uglob + for i in range(ntot_u): + xglob = M_data[p, 1] + yglob = M_data[p, 2] + uglob = M_data[p, 4] + + # Compute analytical value in local coordinate system: + XGLOB = np.array([xglob, yglob]) - displ + XLOC = TROTMAT @ XGLOB - dispxy + + # Analytical Velocity in local axes (u, v): + ULOC = np.array([ + ann_u(XLOC[0], XLOC[1], T), + ann_v(XLOC[0], XLOC[1], T) + ]) + + # Analytical Velocity in global axes: + UGLOB = ROTMAT @ ULOC + uglob_ann = UGLOB[0] + + e_u_vec[i] = uglob - uglob_ann + p += 1 + + # Calculate norms for U + error_lists['e_u_1'].append(np.linalg.norm(e_u_vec, 1) / ntot_u) # L1/N + error_lists['e_u_2'].append(np.linalg.norm(e_u_vec, 2) / np.sqrt(ntot_u)) # L2/sqrt(N) (RMS-like) + error_lists['e_u_i'].append(np.linalg.norm(e_u_vec, np.inf)) # L_inf + + + # --- Second V: Velocity component along y-axis --- + # M_data columns: 0: iscf, 1: xglob, 2: yglob, 3: area, 4: vglob + for i in range(ntot_v): + xglob = M_data[p, 1] + yglob = M_data[p, 2] + vglob = M_data[p, 4] + + # Compute analytical value in local coordinate system: + XGLOB = np.array([xglob, yglob]) - displ + XLOC = TROTMAT @ XGLOB - dispxy + + # Analytical Velocity in local axes (u, v): + ULOC = np.array([ + ann_u(XLOC[0], XLOC[1], T), + ann_v(XLOC[0], XLOC[1], T) + ]) + + # Analytical Velocity in global axes: + UGLOB = ROTMAT @ ULOC + vglob_ann = UGLOB[1] + + e_v_vec[i] = vglob - vglob_ann + p += 1 + + # Calculate norms for V + error_lists['e_v_1'].append(np.linalg.norm(e_v_vec, 1) / ntot_v) + error_lists['e_v_2'].append(np.linalg.norm(e_v_vec, 2) / np.sqrt(ntot_v)) + error_lists['e_v_i'].append(np.linalg.norm(e_v_vec, np.inf)) + + # --- Finally cell centered variables (Z, H, P) --- + # M_data columns: 0: iscf, 1: xglob, 2: yglob, 3: vol, 4: z, 5: H, 6: pres + DELP = 0.0 + p_c = p # Start index for cell-centered variables + + for i in range(ntot_c): + # Data from the CSV row: + xglob = M_data[p, 1] + yglob = M_data[p, 2] + z = M_data[p, 4] + H = M_data[p, 5] + pres = M_data[p, 6] + + # Compute analytical scalars in local coordinate system: + XGLOB = np.array([xglob, yglob]) - displ + XLOC = TROTMAT @ XGLOB - dispxy + + z_ann = ann_z(XLOC[0], XLOC[1], T) + H_ann = ann_H(XLOC[0], XLOC[1], T) + p_ann = ann_p(XLOC[0], XLOC[1], T) + + if i == 0: + # Pressure constant of integration (shift) + DELP = pres - p_ann + + # Error calculation: + e_z_vec[i] = z - z_ann + e_H_vec[i] = H - H_ann - DELP / rho + e_p_vec[i] = pres - p_ann - DELP # Pressure error is corrected by the shift DELP + + p += 1 + + # Calculate norms for Z + error_lists['e_z_1'].append(np.linalg.norm(e_z_vec, 1) / ntot_c) + error_lists['e_z_2'].append(np.linalg.norm(e_z_vec, 2) / np.sqrt(ntot_c)) + error_lists['e_z_i'].append(np.linalg.norm(e_z_vec, np.inf)) + + # Calculate norms for H + error_lists['e_H_1'].append(np.linalg.norm(e_H_vec, 1) / ntot_c) + error_lists['e_H_2'].append(np.linalg.norm(e_H_vec, 2) / np.sqrt(ntot_c)) + error_lists['e_H_i'].append(np.linalg.norm(e_H_vec, np.inf)) + + # Calculate norms for P + error_lists['e_p_1'].append(np.linalg.norm(e_p_vec, 1) / ntot_c) + error_lists['e_p_2'].append(np.linalg.norm(e_p_vec, 2) / np.sqrt(ntot_c)) + error_lists['e_p_i'].append(np.linalg.norm(e_p_vec, np.inf)) + + # Convert error lists to numpy arrays and store in the file entry + for key, val_list in error_lists.items(): + file_entry['errors'][key] = np.array(val_list) + + files[ifile] = file_entry # Update the entry in the list + + +# --- 5. Generate Plots and Warnings --- +print("\n--- Summary of Calculated Errors and Convergence Warnings ---") + +for ifile in range(len(files)): + file_entry = files[ifile] + dx_values = np.array(file_entry['dx']) + + # Axis limits based on MATLAB original + if ifile == 0 or ifile == 1: + axisval=[1e-2, 1, 1e-7, 1e-1] + elif ifile == 2: + axisval=[1e-2, 1, 1e-6, 1e-1] + else: + axisval=[1e-2, 1, 1e-7, 1e-1] + + p_L2_u = np.log(file_entry['errors']['e_u_2'][:-1] / file_entry['errors']['e_u_2'][1:]) / np.log(dx_values[:-1] / dx_values[1:]) + print(f"\nConvergence Order (p) for L2 Error u, {file_entry['nameout']}: {p_L2_u}") + p_L2_z = np.log(file_entry['errors']['e_z_2'][:-1] / file_entry['errors']['e_z_2'][1:]) / np.log(dx_values[:-1] / dx_values[1:]) + print(f"Convergence Order (p) for L2 Error z, {file_entry['nameout']}: {p_L2_z}") + + fig = fdsplotlib.plot_to_fig(x_data=dx_values, y_data=file_entry['errors']['e_z_2'], marker_style='rs-',plot_type='loglog', + revision_label=version_string,x_min=axisval[0],x_max=axisval[1],y_min=axisval[2],y_max=axisval[3], + data_label='FDS $Z$', + x_label='$\\Delta x$ (m)', + y_label='$L_2$ Error') + + fdsplotlib.plot_to_fig(x_data=dx_values, y_data= file_entry['errors']['e_u_2'], marker_style='g>-', + figure_handle=fig, + data_label='FDS $u$') + + fdsplotlib.plot_to_fig(x_data=dx_values, y_data= file_entry['errors']['e_H_2'], marker_style='b+>-', + figure_handle=fig, + data_label='FDS $H$') + + O1_ref = 10**-1 * dx_values + fdsplotlib.plot_to_fig(x_data=dx_values, y_data= O1_ref, marker_style='k--', + figure_handle=fig, + data_label='$O(\\Delta x)$') + O2_ref = 2 * 10**-3 * dx_values**2 + fdsplotlib.plot_to_fig(x_data=dx_values, y_data= O2_ref, marker_style='k-', + figure_handle=fig, + data_label='$O(\\Delta x^2)$') + + plotname = plotdir + file_entry['nameout'] + '.pdf' + plt.savefig(plotname, format='pdf') + plt.close() + + # --- Warning check logic --- + last_z_error = file_entry['errors']['e_z_2'][-1] + last_u_error = file_entry['errors']['e_u_2'][-1] + last_H_error = file_entry['errors']['e_H_2'][-1] + + if ifile == 1: # OBST case (ifile == 2 in MATLAB) + print(f"\n*** Tolerance Check for OBST Case ({file_entry['nameout'].replace('_',' ')}) ***") + if last_z_error > 2e-6: + print(f" Warning (Z): {last_z_error:.2e} > 2e-6 tolerance.") + if last_u_error > 2e-5: + print(f" Warning (u): {last_u_error:.2e} > 2e-5 tolerance.") + if last_H_error > 1.5e-3: + print(f" Warning (H): {last_H_error:.2e} > 1.5e-3 tolerance.") + elif ifile == 2: # 45deg stm case (ifile == 3 in MATLAB) + print(f"\n*** Tolerance Check for 45deg STM Case ({file_entry['nameout'].replace('_',' ')}) ***") + if last_z_error > 6.5e-5: + print(f" Warning (Z): {last_z_error:.2e} > 6.5e-5 tolerance.") + if last_u_error > 8.5e-4: + print(f" Warning (u): {last_u_error:.2e} > 8.5e-4 tolerance.") + if last_H_error > 2.5e-3: + print(f" Warning (H): {last_H_error:.2e} > 2.5e-3 tolerance.") + + diff --git a/Verification/Species/EDC_methane.py b/Verification/Species/EDC_methane.py new file mode 100644 index 00000000000..977e432159e --- /dev/null +++ b/Verification/Species/EDC_methane.py @@ -0,0 +1,385 @@ + +# Generate exact solution for EDC_methane verification cases + +import numpy as np +from scipy.integrate import ode +import csv + +def EDC_spec_methane(t, y): + """ + Calculate derivatives for methane combustion system. + + Species indices: + (0) - nitrogen + (1) - methane + (2) - oxygen + (3) - co2 + (4) - h2o + + Args: + t: Time variable + y: Array of species concentrations + + Returns: + dy: Array of derivatives for each species + """ + + m = 1 # chemical reaction coefficient fuel + s = 2 # chemical reaction coefficient oxygen + c = 1 # chemical reaction coefficient carbon dioxide + h = 2 # chemical reaction coefficient water vapor + tau_mix = 0.125 # set mixing time [s] + + # ----------------------- + # Logic for setting molar production rate + # No suppression in this case + # ----------------------- + r1 = min([y[1], y[2]/s]) + r_co2 = np.zeros(1) + r_co2[0] = m * r1 / tau_mix + + dy = np.zeros(5) + dy[0] = 0.0 + dy[1] = -(1/m) * r_co2[0] + dy[2] = -(s/m) * r_co2[0] + dy[3] = c * r_co2[0] + dy[4] = h * r_co2[0] + + return dy + + +temp_0 = 298.15 # initial temperature [K] +volume = 0.001 # volume [m^3] +rho = 1.1294838 # density [kg/m^3] +R = 8.3145 # gas constant [J/mol K] +mass = rho * volume # mass [kg] + +yf0 = np.array([0.7247698, 0.0551642, 0.220066, 0.0, 0.0]) # initial mass fraction + +y_MW = np.array([28.0134, 16.042460, 31.9988, 44.0095, 18.01528]) # [g/mol] + +y_hf = np.array([0.0, -74873, 0.0, -393522, -241826]) # [J/mol] + +N0 = (1000 * volume * rho * yf0) / y_MW # initial moles + +pres_0 = temp_0 * np.sum(N0) * R / volume # initial pressure + +tspan = np.arange(0, 2.1, 0.1) + +y0 = np.array([N0[0], N0[1], N0[2], N0[3], N0[4]]) + +solver = ode(EDC_spec_methane) +solver.set_integrator('dopri5', rtol=1e-10, atol=1e-10) +solver.set_initial_value(y0, tspan[0]) + +# Integrate over time +T = [tspan[0]] +Y = [y0] +for t in tspan[1:]: + solver.integrate(t) + T.append(solver.t) + Y.append(solver.y) + +T = np.array(T) +Y = np.array(Y) + +Nf = Y[-1, :] + +# Convert back to mass fraction +yff = np.zeros((len(T), 5)) +yff[:, 0] = (y_MW[0] * Y[:, 0]) / (1000 * rho * volume) +yff[:, 1] = (y_MW[1] * Y[:, 1]) / (1000 * rho * volume) +yff[:, 2] = (y_MW[2] * Y[:, 2]) / (1000 * rho * volume) +yff[:, 3] = (y_MW[3] * Y[:, 3]) / (1000 * rho * volume) +yff[:, 4] = (y_MW[4] * Y[:, 4]) / (1000 * rho * volume) + +# Determine final temperature and pressure +tol0 = 1e5 +tol = 0.1 +Tf_guess = np.array([1000.0, 2000.0, 3000.0]) # final temperature guess [K] +count = 0 + +while abs(tol0) > tol: + # --------- + # NOTE: coefficients were found from NIST Webbook + # --------- + + # initial temperature coeffs + coeff0 = np.zeros((7, 5)) + coeff0[:, 0] = [28.98641, 1.853978, -9.647459, 16.63537, 0.000117, -8.671914, 0.0] + coeff0[:, 1] = [-0.703029, 108.4773, -42.52157, 5.862788, 0.678565, -76.84376, -74.87310] + coeff0[:, 2] = [31.32234, -20.23531, 57.86644, -36.50624, -0.007374, -8.903471, 0.0] + coeff0[:, 3] = [24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, -393.5224] + coeff0[:, 4] = [30.09200, 6.832514, 6.793435, -2.534480, 0.082139, -250.8810, -241.8264] + + coeff = np.zeros((7, 5)) + + # nitrogen cp coeffs [J/mol K] + if Tf_guess[1] <= 500: + coeff[:, 0] = [28.98641, 1.853978, -9.647459, 16.63537, 0.000117, -8.671914, 0.0] + elif 500 < Tf_guess[1] <= 2000: + coeff[:, 0] = [19.50583, 19.88705, -8.598535, 1.369784, 0.527601, -4.935202, 0.0] + else: + coeff[:, 0] = [35.51872, 1.128728, -0.196103, 0.014662, -4.553760, -18.97091, 0.0] + + # methane cp coeffs [J/mol K] + if Tf_guess[1] <= 1300: + coeff[:, 1] = [-0.703029, 108.4773, -42.52157, 5.862788, 0.678565, -76.84376, -74.87310] + else: + coeff[:, 1] = [85.81217, 11.26467, -2.114146, 0.138190, -26.42221, -153.5327, -74.87310] + + # oxygen cp coeffs [J/mol K] + if Tf_guess[1] <= 700: + coeff[:, 2] = [31.32234, -20.23531, 57.86644, -36.50624, -0.007374, -8.903471, 0.0] + elif 700 < Tf_guess[1] <= 2000: + coeff[:, 2] = [30.03235, 8.772972, -3.988133, 0.788313, -0.741599, -11.32468, 0.0] + else: + coeff[:, 2] = [20.91111, 10.72071, -2.020498, 0.146449, 9.245722, 5.337651, 0.0] + + # carbon dioxide cp coeffs [J/mol K] + if Tf_guess[1] <= 1200: + coeff[:, 3] = [24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, -393.5224] + else: + coeff[:, 3] = [58.16639, 2.720074, -0.492289, 0.038844, -6.447293, -425.9186, -393.5224] + + # water vapor cp coeffs [J/mol K] + if Tf_guess[1] <= 1700: + coeff[:, 4] = [30.09200, 6.832514, 6.793435, -2.534480, 0.082139, -250.8810, -241.8264] + else: + coeff[:, 4] = [41.96246, 8.622053, -1.499780, 0.098199, -11.15764, -272.1797, -241.8264] + + t = Tf_guess / 1000 + t0 = temp_0 / 1000 + del_h = np.zeros((5, 3)) + + for j in range(3): + for i in range(5): + del_h[i, j] = ((coeff[0, i] * t[j] + (1/2) * coeff[1, i] * t[j]**2 + (1/3) * coeff[2, i] * t[j]**3 + + (1/4) * coeff[3, i] * t[j]**4 - coeff[4, i] / t[j] + coeff[5, i] - coeff[6, i]) - + (coeff0[0, i] * t0 + (1/2) * coeff0[1, i] * t0**2 + (1/3) * coeff0[2, i] * t0**3 + + (1/4) * coeff0[3, i] * t0**4 - coeff0[4, i] / t0 + coeff0[5, i] - coeff0[6, i])) * 1000 + + h_fc = Nf[3] * y_hf[3] + Nf[4] * y_hf[4] - N0[1] * y_hf[1] + del_hN = Nf[0] * del_h[0, :] + Nf[3] * del_h[3, :] + Nf[4] * del_h[4, :] + RT0 = np.sum(N0) * R * temp_0 + RTf = -np.sum(Nf) * R * Tf_guess + + Tf_solve = h_fc + del_hN + RT0 + RTf + tol0 = Tf_guess[0] - Tf_guess[2] + + if Tf_solve[1] < 0: + if Tf_solve[0] > 0: + mp = 0.5 * (Tf_guess[0] + Tf_guess[1]) + Tf_guess = np.array([Tf_guess[0], mp, Tf_guess[1]]) + else: + mp = 0.5 * (Tf_guess[1] + Tf_guess[2]) + Tf_guess = np.array([Tf_guess[1], mp, Tf_guess[2]]) + else: + if Tf_solve[0] < 0: + mp = 0.5 * (Tf_guess[0] + Tf_guess[1]) + Tf_guess = np.array([Tf_guess[0], mp, Tf_guess[1]]) + else: + mp = 0.5 * (Tf_guess[1] + Tf_guess[2]) + Tf_guess = np.array([Tf_guess[1], mp, Tf_guess[2]]) + + count = count + 1 + +Tf = Tf_guess[1] +Pf = (Tf_guess[1] * np.sum(Nf) * R) / volume +dP = Pf - pres_0 +Tf = Tf - 273.15 + +tss = np.array([1.0, 1.25, 1.5, 1.75, 2.0]) +Tf_vec = np.array([Tf, Tf, Tf, Tf, Tf]) +dP_vec = np.array([dP, dP, dP, dP, dP]) + +yff_output = np.column_stack((tspan, yff[:, 2], yff[:, 1], yff[:, 3], yff[:, 4])) + +header1 = ['Time', 'CH4', 'O2', 'CO2', 'H2O'] +with open('reactionrate_EDC_flim_1step_CH4_spec.csv', 'w', newline='') as fid: + writer = csv.writer(fid) + writer.writerow(header1) + for j in range(len(tspan)): + formatted_row = [f"{value:.6f}" for value in yff_output[j, :]] + writer.writerow(formatted_row) + +header1 = ['Time', 'TEMP', 'PRES'] +with open('reactionrate_EDC_flim_1step_CH4_temppres.csv', 'w', newline='') as fid: + writer = csv.writer(fid) + writer.writerow(header1) + for j in range(len(tss)): + writer.writerow([f"{tss[j]:.6f}", f"{Tf_vec[j]:.6f}", f"{dP_vec[j]:.6f}"]) + + + + +# Non premix case/fuel rich + +temp_0 = 298.15 # initial temperature [K] +volume = 0.001 # volume [m^3] +rho = 1.1294838 # density [kg/m^3] +R = 8.3145 # gas constant [J/mol K] +mass = rho * volume # mass [kg] + +yf0 = np.array([0.7252, 0.0582, 0.2166, 0.0, 0.0]) +y_MW = np.array([28.0134, 16.042460, 31.9988, 44.0095, 18.01528]) # [g/mol] + +y_hf = np.array([0.0, -74873, 0.0, -393522, -241826]) # [J/mol] + +N0 = (1000 * volume * rho * yf0) / y_MW # initial moles + +pres_0 = temp_0 * np.sum(N0) * R / volume # initial pressure + +tspan = np.arange(0, 65, 5) + +y0 = np.array([N0[0], N0[1], N0[2], N0[3], N0[4]]) + +solver = ode(EDC_spec_methane) +#solver.set_integrator('dopri5', rtol=1e-10, atol=1e-10) +solver.set_integrator('vode', method='bdf') +solver.set_initial_value(y0, tspan[0]) + +Y = [] +T = [] +for t_val in tspan[1:]: + if solver.successful(): + solver.integrate(t_val) + Y.append(solver.y.copy()) + T.append(t_val) + else: + break + +Y = np.array(Y) +T = np.array(T) + +Nf = Y[-1, :] + +# Convert back to mass fraction +yff = np.zeros((len(tspan), 5)) +yff[1:, 0] = (y_MW[0] * Y[:, 0]) / (1000 * rho * volume) +yff[1:, 1] = (y_MW[1] * Y[:, 1]) / (1000 * rho * volume) +yff[1:, 2] = (y_MW[2] * Y[:, 2]) / (1000 * rho * volume) +yff[1:, 3] = (y_MW[3] * Y[:, 3]) / (1000 * rho * volume) +yff[1:, 4] = (y_MW[4] * Y[:, 4]) / (1000 * rho * volume) + +# Determine final temperature and pressure +tol0 = 1e5 +tol = 0.1 +Tf_guess = np.array([1000.0, 2000.0, 3000.0]) # final temperature guess [K] +count = 0 + +while abs(tol0) > tol: + # --------- + # NOTE: coefficients were found from NIST Webbook + # --------- + + # initial temperature coeffs + coeff0 = np.zeros((7, 5)) + coeff0[:, 0] = np.array([28.98641, 1.853978, -9.647459, 16.63537, 0.000117, -8.671914, 0.0]) + coeff0[:, 1] = np.array([-0.703029, 108.4773, -42.52157, 5.862788, 0.678565, -76.84376, -74.87310]) + coeff0[:, 2] = np.array([31.32234, -20.23531, 57.86644, -36.50624, -0.007374, -8.903471, 0.0]) + coeff0[:, 3] = np.array([24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, -393.5224]) + coeff0[:, 4] = np.array([30.09200, 6.832514, 6.793435, -2.534480, 0.082139, -250.8810, -241.8264]) + + coeff = np.zeros((7, 5)) + + # nitrogen cp coeffs [J/mol K] + if Tf_guess[1] <= 500: + coeff[:, 0] = np.array([28.98641, 1.853978, -9.647459, 16.63537, 0.000117, -8.671914, 0.0]) + elif 500 < Tf_guess[1] <= 2000: + coeff[:, 0] = np.array([19.50583, 19.88705, -8.598535, 1.369784, 0.527601, -4.935202, 0.0]) + else: + coeff[:, 0] = np.array([35.51872, 1.128728, -0.196103, 0.014662, -4.553760, -18.97091, 0.0]) + + # methane cp coeffs [J/mol K] + if Tf_guess[1] <= 1300: + coeff[:, 1] = np.array([-0.703029, 108.4773, -42.52157, 5.862788, 0.678565, -76.84376, -74.87310]) + else: + coeff[:, 1] = np.array([85.81217, 11.26467, -2.114146, 0.138190, -26.42221, -153.5327, -74.87310]) + + # oxygen cp coeffs [J/mol K] + if Tf_guess[1] <= 700: + coeff[:, 2] = np.array([31.32234, -20.23531, 57.86644, -36.50624, -0.007374, -8.903471, 0.0]) + elif 700 < Tf_guess[1] <= 2000: + coeff[:, 2] = np.array([30.03235, 8.772972, -3.988133, 0.788313, -0.741599, -11.32468, 0.0]) + else: + coeff[:, 2] = np.array([20.91111, 10.72071, -2.020498, 0.146449, 9.245722, 5.337651, 0.0]) + + # carbon dioxide cp coeffs [J/mol K] + if Tf_guess[1] <= 1200: + coeff[:, 3] = np.array([24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, -393.5224]) + else: + coeff[:, 3] = np.array([58.16639, 2.720074, -0.492289, 0.038844, -6.447293, -425.9186, -393.5224]) + + # water vapor cp coeffs [J/mol K] + if Tf_guess[1] <= 1700: + coeff[:, 4] = np.array([30.09200, 6.832514, 6.793435, -2.534480, 0.082139, -250.8810, -241.8264]) + else: + coeff[:, 4] = np.array([41.96246, 8.622053, -1.499780, 0.098199, -11.15764, -272.1797, -241.8264]) + + t = Tf_guess / 1000 + t0 = temp_0 / 1000 + + del_h = np.zeros((5, 3)) + for j in range(3): + for i in range(5): + del_h[i, j] = ((coeff[0, i] * t[j] + (1/2) * coeff[1, i] * t[j]**2 + (1/3) * coeff[2, i] * t[j]**3 + + (1/4) * coeff[3, i] * t[j]**4 - coeff[4, i] / t[j] + coeff[5, i] - coeff[6, i]) - + (coeff0[0, i] * t0 + (1/2) * coeff0[1, i] * t0**2 + (1/3) * coeff0[2, i] * t0**3 + + (1/4) * coeff0[3, i] * t0**4 - coeff0[4, i] / t0 + coeff0[5, i] - coeff0[6, i])) * 1000 + + h_fc = Nf[3] * y_hf[3] + Nf[4] * y_hf[4] - N0[1] * y_hf[1] + del_hN = Nf[0] * del_h[0, :] + Nf[3] * del_h[3, :] + Nf[4] * del_h[4, :] + Nf[1] * del_h[1, :] + RT0 = np.sum(N0) * R * temp_0 + RTf = -np.sum(Nf) * R * Tf_guess + + Tf_solve = h_fc + del_hN + RT0 + RTf + tol0 = Tf_guess[0] - Tf_guess[2] + + if Tf_solve[1] < 0: + if Tf_solve[0] > 0: + mp = 0.5 * (Tf_guess[0] + Tf_guess[1]) + Tf_guess = np.array([Tf_guess[0], mp, Tf_guess[1]]) + else: + mp = 0.5 * (Tf_guess[1] + Tf_guess[2]) + Tf_guess = np.array([Tf_guess[1], mp, Tf_guess[2]]) + else: + if Tf_solve[0] < 0: + mp = 0.5 * (Tf_guess[0] + Tf_guess[1]) + Tf_guess = np.array([Tf_guess[0], mp, Tf_guess[1]]) + else: + mp = 0.5 * (Tf_guess[1] + Tf_guess[2]) + Tf_guess = np.array([Tf_guess[1], mp, Tf_guess[2]]) + + count = count + 1 + +Tf = Tf_guess[1] +Pf = (Tf_guess[1] * np.sum(Nf) * R) / volume +dP = Pf - pres_0 +Tf = Tf - 273.15 + +tss = 5 * (np.arange(1, 17)) + 25 +Tf_vec = Tf * np.ones(16) +dP_vec = dP * np.ones(16) + +tss = np.array([35, 40, 45, 50, 55, 60]) +Tf_final = np.array([Tf, Tf, Tf, Tf, Tf, Tf]) +dP_final = np.array([dP, dP, dP, dP, dP, dP]) + +yff[:, 0] = tspan + +header1 = ['Time', 'CH4', 'O2', 'CO2', 'H2O'] +with open('reactionrate_EDC_1step_CH4_nonmix_spec.csv', 'w', newline='') as fid: + writer = csv.writer(fid) + writer.writerow(header1) + for j in range(7, len(tspan)): + writer.writerow([f"{yff[j,0]:.6f}", f"{yff[j,1]:.6f}", f"{yff[j,2]:.6f}", f"{yff[j,3]:.6f}", f"{yff[j, 4]:.6f}"]) + +header1 = ['Time', 'TEMP', 'PRES'] +with open('reactionrate_EDC_1step_CH4_nonmix_temppres.csv', 'w', newline='') as fid: + writer = csv.writer(fid) + writer.writerow(header1) + for j in range(len(tss)): + writer.writerow([f"{tss[j]:.6f}", f"{Tf_final[j]:.6f}", f"{dP_final[j]:.6f}"]) + + diff --git a/Verification/Species/EDC_propane.py b/Verification/Species/EDC_propane.py new file mode 100644 index 00000000000..2dcbfa607d7 --- /dev/null +++ b/Verification/Species/EDC_propane.py @@ -0,0 +1,209 @@ + +# Verification of EDM for propane + +import numpy as np +from scipy.integrate import ode +import csv + +def EDC_spec_propane(t, y): + """ + ODE system for propane combustion with EDC (Eddy Dissipation Concept) + + Species: + (0) - nitrogen + (1) - propane + (2) - oxygen + (3) - co2 + (4) - h2o + """ + + m = 1 # chemical reaction coefficient fuel + s = 5 # chemical reaction coefficient oxygen + c = 3 # chemical reaction coefficient carbon dioxide + h = 4 # chemical reaction coefficient water vapor + tau_mix = 0.125 # set mixing time [s] + + r1 = min([y[1], y[2]/s]) + r_co2 = np.zeros(1) + r_co2[0] = m * r1 / tau_mix + + dy = np.zeros(5) + dy[0] = 0.0 + dy[1] = -(1/m) * r_co2[0] + dy[2] = -(s/m) * r_co2[0] + dy[3] = c * r_co2[0] + dy[4] = h * r_co2[0] + + return dy + + +temp_0 = 298.15 # initial temperature [K] +volume = 0.001 # volume [m^3] +rho = 1.2043424 # density [kg/m^3] +R = 8.3145 # gas constant [J/mol K] +mass = rho * volume # mass [kg] + +yf0 = np.array([0.720828, 0.060321, 0.218851, 0.0, 0.0]) # initial mass fraction +y_MW = np.array([28.0134, 44.09562, 31.9988, 44.0095, 18.01528]) # [g/mol] +y_hf = np.array([0.0, -103850, 0.0, -393522, -241826]) # [J/mol] + +N0 = (1000 * volume * rho * yf0) / y_MW # initial moles + +pres_0 = temp_0 * np.sum(N0) * R / volume # initial pressure + +tspan = np.arange(0, 2.1, 0.1) +y0 = np.array([N0[0], N0[1], N0[2], N0[3], N0[4]]) + +# Setup ODE solver with tight tolerances +solver = ode(EDC_spec_propane) +solver.set_integrator('vode', method='adams', rtol=1e-10, atol=1e-10, max_step=0.1) +solver.set_initial_value(y0, tspan[0]) + +# Integrate over time span +Y = [] +T = [] +for t in tspan[1:]: + if solver.successful(): + solver.integrate(t) + Y.append(solver.y.copy()) + T.append(t) + else: + break + +Y = np.array(Y) +T = np.array(T) + +Nf = Y[-1, :] + +# Convert back to mass fraction +yff = np.zeros((len(T), 6)) +yff[:, 0] = T +yff[:, 1] = (y_MW[2] * Y[:, 2]) / (1000 * rho * volume) # O2 +yff[:, 2] = (y_MW[1] * Y[:, 1]) / (1000 * rho * volume) # C3H8 +yff[:, 3] = (y_MW[3] * Y[:, 3]) / (1000 * rho * volume) # CO2 +yff[:, 4] = (y_MW[4] * Y[:, 4]) / (1000 * rho * volume) # H2O + +# ----------------------- +# Determine final temperature and pressure +# ----------------------- +tol0 = 1e5 +tol = 0.1 +Tf_guess = np.array([1000.0, 2000.0, 3000.0]) # final temperature guess [K] +count = 0 + +while abs(tol0) > tol: + # --------- + # NOTE: coefficients were found from NIST Webbook + # --------- + + # initial temperature coeffs + coeff0 = np.zeros((7, 5)) + coeff0[:, 0] = np.array([28.98641, 1.853978, -9.647459, 16.63537, 0.000117, -8.671914, 0.0]) + coeff0[:, 1] = np.array([1, 1, 1, 1, 1, 1, 1]) # placeholder for propane + coeff0[:, 2] = np.array([31.32234, -20.23531, 57.86644, -36.50624, -0.007374, -8.903471, 0.0]) + coeff0[:, 3] = np.array([24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, -393.5224]) + coeff0[:, 4] = np.array([30.09200, 6.832514, 6.793435, -2.534480, 0.082139, -250.8810, -241.8264]) + + # Initialize coefficient array + coeff = np.zeros((7, 5)) + + # nitrogen cp coeffs [J/mol K] + if Tf_guess[1] <= 500: + coeff[:, 0] = np.array([28.98641, 1.853978, -9.647459, 16.63537, 0.000117, -8.671914, 0.0]) + elif 500 < Tf_guess[1] <= 2000: + coeff[:, 0] = np.array([19.50583, 19.88705, -8.598535, 1.369784, 0.527601, -4.935202, 0.0]) + else: + coeff[:, 0] = np.array([35.51872, 1.128728, -0.196103, 0.014662, -4.553760, -18.97091, 0.0]) + + # propane cp coeffs [J/mol K] + if Tf_guess[1] <= 1300: + coeff[:, 1] = np.array([1, 1, 1, 1, 1, 1, 1]) # placeholder for propane + else: + coeff[:, 1] = np.array([1, 1, 1, 1, 1, 1, 1]) # placeholder for propane + + # oxygen cp coeffs [J/mol K] + if Tf_guess[1] <= 700: + coeff[:, 2] = np.array([31.32234, -20.23531, 57.86644, -36.50624, -0.007374, -8.903471, 0.0]) + elif 700 < Tf_guess[1] <= 2000: + coeff[:, 2] = np.array([30.03235, 8.772972, -3.988133, 0.788313, -0.741599, -11.32468, 0.0]) + else: + coeff[:, 2] = np.array([20.91111, 10.72071, -2.020498, 0.146449, 9.245722, 5.337651, 0.0]) + + # carbon dioxide cp coeffs [J/mol K] + if Tf_guess[1] <= 1200: + coeff[:, 3] = np.array([24.99735, 55.18696, -33.69137, 7.948387, -0.136638, -403.6075, -393.5224]) + else: + coeff[:, 3] = np.array([58.16639, 2.720074, -0.492289, 0.038844, -6.447293, -425.9186, -393.5224]) + + # water vapor cp coeffs [J/mol K] + if Tf_guess[1] <= 1700: + coeff[:, 4] = np.array([30.09200, 6.832514, 6.793435, -2.534480, 0.082139, -250.8810, -241.8264]) + else: + coeff[:, 4] = np.array([41.96246, 8.622053, -1.499780, 0.098199, -11.15764, -272.1797, -241.8264]) + + t = Tf_guess / 1000 + t0 = temp_0 / 1000 + + del_h = np.zeros((5, 3)) + for j in range(3): + for i in range(5): + del_h[i, j] = ((coeff[0, i] * t[j] + (1/2) * coeff[1, i] * t[j]**2 + (1/3) * coeff[2, i] * t[j]**3 + (1/4) * coeff[3, i] * t[j]**4 - coeff[4, i] / t[j] + coeff[5, i] - coeff[6, i]) \ + - (coeff0[0, i] * t0 + (1/2) * coeff0[1, i] * t0**2 + (1/3) * coeff0[2, i] * t0**3 + (1/4) * coeff0[3, i] * t0**4 - coeff0[4, i] / t0 + coeff0[5, i] - coeff0[6, i])) \ + * 1000 + + # overwrite del_h[1, :] for propane (index 1 corresponds to propane) + coeff2 = np.array([1.67536e-9, -5.46675e-6, 4.38029e-3, 2.7949, 6.11728e2]) + for j in range(3): + del_h[1, j] = ((1/5) * coeff2[0] * Tf_guess[j]**5 + (1/4) * coeff2[1] * Tf_guess[j]**4 + (1/3) * coeff2[2] * Tf_guess[j]**3 + (1/2) * coeff2[3] * Tf_guess[j]**2 + coeff2[4] * Tf_guess[j]) \ + - ((1/5) * coeff2[0] * temp_0**5 + (1/4) * coeff2[1] * temp_0**4 + (1/3) * coeff2[2] * temp_0**3 + (1/2) * coeff2[3] * temp_0**2 + coeff2[4] * temp_0) \ + * (Tf_guess[j] - temp_0) + + h_fc = Nf[3] * y_hf[3] + Nf[4] * y_hf[4] - N0[1] * y_hf[1] + del_hN = Nf[0] * del_h[0, :] + Nf[3] * del_h[3, :] + Nf[4] * del_h[4, :] + RT0 = np.sum(N0) * R * temp_0 + RTf = -np.sum(Nf) * R * Tf_guess + + Tf_solve = h_fc + del_hN + RT0 + RTf + tol0 = Tf_guess[0] - Tf_guess[2] + + if Tf_solve[1] < 0: + if Tf_solve[0] > 0: + mp = 0.5 * (Tf_guess[0] + Tf_guess[1]) + Tf_guess = np.array([Tf_guess[0], mp, Tf_guess[1]]) + else: + mp = 0.5 * (Tf_guess[1] + Tf_guess[2]) + Tf_guess = np.array([Tf_guess[1], mp, Tf_guess[2]]) + else: + if Tf_solve[0] < 0: + mp = 0.5 * (Tf_guess[0] + Tf_guess[1]) + Tf_guess = np.array([Tf_guess[0], mp, Tf_guess[1]]) + else: + mp = 0.5 * (Tf_guess[1] + Tf_guess[2]) + Tf_guess = np.array([Tf_guess[1], mp, Tf_guess[2]]) + + count = count + 1 + +Tf = Tf_guess[1] +Pf = (Tf_guess[1] * np.sum(Nf) * R) / volume +dP = Pf - pres_0 +Tf = Tf - 273.15 + +tss = np.array([1, 1.25, 1.5, 1.75, 2]) +Tf_vec = np.array([Tf, Tf, Tf, Tf, Tf]) +dP_vec = np.array([dP, dP, dP, dP, dP]) + +header1 = ['Time', 'O2', 'C3H8', 'CO2', 'H2O'] +filename1 = '../../Verification/Species/reactionrate_EDC_flim_1step_C3H8_spec.csv' +with open('reactionrate_EDC_flim_1step_C3H8_spec.csv', 'w', newline='') as fid: + writer = csv.writer(fid) + writer.writerow(header1) + for j in range(len(T)): + writer.writerow([f"{yff[j,0]:.6f}", f"{yff[j,1]:.6f}", f"{yff[j,2]:.6f}", f"{yff[j,3]:.6f}", f"{yff[j,4]:.6f}"]) + +header1 = ['Time', 'TEMP', 'PRES'] +with open('reactionrate_EDC_flim_1step_C3H8_temppres.csv', 'w', newline='') as fid: + writer = csv.writer(fid) + writer.writerow(header1) + for j in range(len(tss)): + writer.writerow([f"{tss[j]:.6f}", f"{Tf_vec[j]:.6f}", f"{dP_vec[j]:.6f}"]) + diff --git a/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_spec.csv b/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_spec.csv index 6eb1775a73c..692f6a298b0 100644 --- a/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_spec.csv +++ b/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_spec.csv @@ -1,7 +1,7 @@ -Time, CH4, O2, CO2, H2O -35.000000, 0.003904, 0.000000, 0.148950, 0.121945 -40.000000, 0.003904, 0.000000, 0.148950, 0.121945 -45.000000, 0.003904, 0.000000, 0.148950, 0.121945 -50.000000, 0.003904, 0.000000, 0.148950, 0.121945 -55.000000, 0.003904, 0.000000, 0.148950, 0.121946 -60.000000, 0.003904, 0.000000, 0.148950, 0.121945 +Time,CH4,O2,CO2,H2O +35.000000,0.003904,-0.000000,0.148950,0.121945 +40.000000,0.003904,0.000000,0.148950,0.121945 +45.000000,0.003904,0.000000,0.148950,0.121945 +50.000000,0.003904,-0.000000,0.148950,0.121945 +55.000000,0.003904,0.000000,0.148950,0.121945 +60.000000,0.003904,0.000000,0.148950,0.121945 diff --git a/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_temppres.csv b/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_temppres.csv index 6079102456e..71f061f85b9 100644 --- a/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_temppres.csv +++ b/Verification/Species/reactionrate_EDC_1step_CH4_nonmix_temppres.csv @@ -1,7 +1,7 @@ -Time, TEMP, PRES -35.000000, 2472.714868, 834061.280941 -40.000000, 2472.714868, 834061.280941 -45.000000, 2472.714868, 834061.280941 -50.000000, 2472.714868, 834061.280941 -55.000000, 2472.714868, 834061.280941 -60.000000, 2472.714868, 834061.280941 +Time,TEMP,PRES +35.000000,2476.956812,835506.727369 +40.000000,2476.956812,835506.727369 +45.000000,2476.956812,835506.727369 +50.000000,2476.956812,835506.727369 +55.000000,2476.956812,835506.727369 +60.000000,2476.956812,835506.727369 diff --git a/Verification/Species/reactionrate_EDC_flim_1step_C3H8_spec.csv b/Verification/Species/reactionrate_EDC_flim_1step_C3H8_spec.csv index 146553045d5..7d641f40ad4 100644 --- a/Verification/Species/reactionrate_EDC_flim_1step_C3H8_spec.csv +++ b/Verification/Species/reactionrate_EDC_flim_1step_C3H8_spec.csv @@ -1,22 +1,21 @@ -Time, O2, C3H8, CO2, H2O -0.000000, 0.060321, 0.218851, 0.000000, 0.000000 -0.100000, 0.027106, 0.098336, 0.099450, 0.054280 -0.200000, 0.012182, 0.044185, 0.144136, 0.078669 -0.300000, 0.005476, 0.019854, 0.164214, 0.089628 -0.400000, 0.002463, 0.008921, 0.173236, 0.094552 -0.500000, 0.001109, 0.004008, 0.177290, 0.096765 -0.600000, 0.000500, 0.001801, 0.179112, 0.097759 -0.700000, 0.000227, 0.000809, 0.179930, 0.098206 -0.800000, 0.000104, 0.000364, 0.180298, 0.098406 -0.900000, 0.000049, 0.000163, 0.180463, 0.098497 -1.000000, 0.000024, 0.000073, 0.180537, 0.098537 -1.100000, 0.000013, 0.000033, 0.180571, 0.098555 -1.200000, 0.000008, 0.000015, 0.180586, 0.098564 -1.300000, 0.000006, 0.000007, 0.180592, 0.098567 -1.400000, 0.000005, 0.000003, 0.180595, 0.098569 -1.500000, 0.000004, 0.000001, 0.180597, 0.098570 -1.600000, 0.000004, 0.000001, 0.180597, 0.098570 -1.700000, 0.000004, 0.000000, 0.180598, 0.098570 -1.800000, 0.000004, 0.000000, 0.180598, 0.098570 -1.900000, 0.000004, 0.000000, 0.180598, 0.098570 -2.000000, 0.000004, 0.000000, 0.180598, 0.098570 +Time,O2,C3H8,CO2,H2O +0.100000,0.098336,0.027106,0.099450,0.054280 +0.200000,0.044185,0.012182,0.144136,0.078669 +0.300000,0.019854,0.005476,0.164214,0.089628 +0.400000,0.008921,0.002463,0.173236,0.094552 +0.500000,0.004008,0.001109,0.177290,0.096765 +0.600000,0.001801,0.000500,0.179112,0.097759 +0.700000,0.000809,0.000227,0.179930,0.098206 +0.800000,0.000364,0.000104,0.180298,0.098406 +0.900000,0.000163,0.000049,0.180463,0.098497 +1.000000,0.000073,0.000024,0.180537,0.098537 +1.100000,0.000033,0.000013,0.180571,0.098555 +1.200000,0.000015,0.000008,0.180586,0.098564 +1.300000,0.000007,0.000006,0.180592,0.098567 +1.400000,0.000003,0.000005,0.180595,0.098569 +1.500000,0.000001,0.000004,0.180597,0.098570 +1.600000,0.000001,0.000004,0.180597,0.098570 +1.700000,0.000000,0.000004,0.180598,0.098570 +1.800000,0.000000,0.000004,0.180598,0.098570 +1.900000,0.000000,0.000004,0.180598,0.098570 +2.000000,0.000000,0.000004,0.180598,0.098570 diff --git a/Verification/Species/reactionrate_EDC_flim_1step_C3H8_temppres.csv b/Verification/Species/reactionrate_EDC_flim_1step_C3H8_temppres.csv index ac2ff3b5dff..44de703832a 100644 --- a/Verification/Species/reactionrate_EDC_flim_1step_C3H8_temppres.csv +++ b/Verification/Species/reactionrate_EDC_flim_1step_C3H8_temppres.csv @@ -1,6 +1,6 @@ -Time, TEMP, PRES -1.000000, 2626.248804, 923739.801851 -1.250000, 2626.248804, 923739.801851 -1.500000, 2626.248804, 923739.801851 -1.750000, 2626.248804, 923739.801851 -2.000000, 2626.248804, 923739.801851 +Time,TEMP,PRES +1.000000,2634.213892,926555.811283 +1.250000,2634.213892,926555.811283 +1.500000,2634.213892,926555.811283 +1.750000,2634.213892,926555.811283 +2.000000,2634.213892,926555.811283 diff --git a/Verification/Species/reactionrate_EDC_flim_1step_CH4_spec.csv b/Verification/Species/reactionrate_EDC_flim_1step_CH4_spec.csv index d1fbf6e3e29..715507876ab 100644 --- a/Verification/Species/reactionrate_EDC_flim_1step_CH4_spec.csv +++ b/Verification/Species/reactionrate_EDC_flim_1step_CH4_spec.csv @@ -1,22 +1,22 @@ -Time, O2, CH4, CO2, H2O -0.000000, 0.055164, 0.220066, 0.000000, 0.000000 -0.100000, 0.024787, 0.098883, 0.083335, 0.068226 -0.200000, 0.011137, 0.044432, 0.120779, 0.098882 -0.300000, 0.005004, 0.019965, 0.137604, 0.112656 -0.400000, 0.002249, 0.008972, 0.145164, 0.118846 -0.500000, 0.001010, 0.004032, 0.148561, 0.121627 -0.600000, 0.000454, 0.001813, 0.150087, 0.122876 -0.700000, 0.000204, 0.000815, 0.150773, 0.123438 -0.800000, 0.000092, 0.000367, 0.151081, 0.123690 -0.900000, 0.000041, 0.000166, 0.151220, 0.123804 -1.000000, 0.000019, 0.000075, 0.151282, 0.123854 -1.100000, 0.000008, 0.000035, 0.151310, 0.123877 -1.200000, 0.000004, 0.000016, 0.151322, 0.123888 -1.300000, 0.000002, 0.000008, 0.151328, 0.123892 -1.400000, 0.000001, 0.000004, 0.151331, 0.123894 -1.500000, 0.000000, 0.000003, 0.151332, 0.123895 -1.600000, 0.000000, 0.000002, 0.151332, 0.123896 -1.700000, 0.000000, 0.000002, 0.151333, 0.123896 -1.800000, 0.000000, 0.000002, 0.151333, 0.123896 -1.900000, 0.000000, 0.000002, 0.151333, 0.123896 -2.000000, 0.000000, 0.000001, 0.151333, 0.123896 +Time,CH4,O2,CO2,H2O +0.000000,0.220066,0.055164,0.000000,0.000000 +0.100000,0.098883,0.024787,0.083335,0.068226 +0.200000,0.044432,0.011137,0.120779,0.098882 +0.300000,0.019965,0.005004,0.137604,0.112656 +0.400000,0.008972,0.002249,0.145164,0.118846 +0.500000,0.004032,0.001010,0.148561,0.121627 +0.600000,0.001813,0.000454,0.150087,0.122876 +0.700000,0.000815,0.000204,0.150773,0.123438 +0.800000,0.000367,0.000092,0.151081,0.123690 +0.900000,0.000166,0.000041,0.151220,0.123804 +1.000000,0.000075,0.000019,0.151282,0.123854 +1.100000,0.000035,0.000008,0.151310,0.123877 +1.200000,0.000016,0.000004,0.151322,0.123888 +1.300000,0.000008,0.000002,0.151328,0.123892 +1.400000,0.000004,0.000001,0.151331,0.123894 +1.500000,0.000003,0.000000,0.151332,0.123895 +1.600000,0.000002,0.000000,0.151332,0.123896 +1.700000,0.000002,0.000000,0.151333,0.123896 +1.800000,0.000002,0.000000,0.151333,0.123896 +1.900000,0.000002,0.000000,0.151333,0.123896 +2.000000,0.000001,0.000000,0.151333,0.123896 diff --git a/Verification/Species/reactionrate_EDC_flim_1step_CH4_temppres.csv b/Verification/Species/reactionrate_EDC_flim_1step_CH4_temppres.csv index 291cdaf8cc3..a6deb4013fd 100644 --- a/Verification/Species/reactionrate_EDC_flim_1step_CH4_temppres.csv +++ b/Verification/Species/reactionrate_EDC_flim_1step_CH4_temppres.csv @@ -1,6 +1,6 @@ -Time, TEMP, PRES -1.000000, 2542.966333, 855722.958695 -1.250000, 2542.966333, 855722.958695 -1.500000, 2542.966333, 855722.958695 -1.750000, 2542.966333, 855722.958695 -2.000000, 2542.966333, 855722.958695 +Time,TEMP,PRES +1.000000,2548.642603,857652.021087 +1.250000,2548.642603,857652.021087 +1.500000,2548.642603,857652.021087 +1.750000,2548.642603,857652.021087 +2.000000,2548.642603,857652.021087 diff --git a/Verification/Sprinklers_and_Sprays/flat_fire.csv b/Verification/Sprinklers_and_Sprays/flat_fire.csv index 8eeeb4e70bc..07492b2a43e 100644 --- a/Verification/Sprinklers_and_Sprays/flat_fire.csv +++ b/Verification/Sprinklers_and_Sprays/flat_fire.csv @@ -1,35 +1,35 @@ Time,x,z,u,w -0.00,0.00000,8.00000,400.00000,0.00000 -0.05,15.06500,7.98970,232.56000,-0.38744 -0.10,24.77800,7.96260,163.93000,-0.69082 -0.15,31.96000,7.92100,126.58000,-0.96759 -0.20,37.66200,7.86600,103.09000,-1.23260 -0.25,42.39000,7.79790,86.95700,-1.49130 -0.30,46.43000,7.71690,75.18800,-1.74630 -0.35,49.95600,7.62330,66.22500,-1.99890 -0.40,53.08400,7.51700,59.17200,-2.24990 -0.45,55.89500,7.39830,53.47600,-2.49980 -0.50,58.44800,7.26710,48.78000,-2.74880 -0.55,60.78600,7.12340,44.84300,-2.99710 -0.60,62.94200,6.96740,41.49400,-3.24500 -0.65,64.94300,6.79890,38.61000,-3.49240 -0.70,66.80900,6.61810,36.10100,-3.73960 -0.75,68.55800,6.42500,33.89800,-3.98640 -0.80,70.20400,6.21950,31.94900,-4.23310 -0.85,71.75700,6.00170,30.21100,-4.47960 -0.90,73.22800,5.77150,28.65300,-4.72590 -0.95,74.62500,5.52910,27.24800,-4.97210 -1.00,75.95500,5.27430,25.97400,-5.21820 -1.05,77.22400,5.00730,24.81400,-5.46420 -1.10,78.43800,4.72790,23.75300,-5.71010 -1.15,79.60100,4.43630,22.77900,-5.95590 -1.20,80.71700,4.13230,21.88200,-6.20170 -1.25,81.79000,3.81610,21.05300,-6.44740 -1.30,82.82300,3.48760,20.28400,-6.69300 -1.35,83.81900,3.14680,19.56900,-6.93860 -1.40,84.78100,2.79370,18.90400,-7.18420 -1.45,85.71000,2.42840,18.28200,-7.42970 -1.50,86.61000,2.05080,17.69900,-7.67520 -1.55,87.48100,1.66090,17.15300,-7.92070 -1.60,88.32600,1.25870,16.63900,-8.16610 -1.65,89.14500,0.84425,16.15500,-8.41150 +0.00000,0.00000,8.00000,400.00000,0.00000 +0.05000,15.06456,7.98968,232.55814,-0.38744 +0.10000,24.77772,7.96255,163.93443,-0.69082 +0.15000,31.96033,7.92102,126.58228,-0.96759 +0.20000,37.66209,7.86598,103.09278,-1.23258 +0.25000,42.39045,7.79787,86.95652,-1.49130 +0.30000,46.42981,7.71691,75.18797,-1.74632 +0.35000,49.95567,7.62327,66.22517,-1.99894 +0.40000,53.08397,7.51705,59.17160,-2.24994 +0.45000,55.89536,7.39830,53.47594,-2.49979 +0.50000,58.44817,7.26708,48.78049,-2.74878 +0.55000,60.78600,7.12343,44.84305,-2.99713 +0.60000,62.94225,6.96738,41.49378,-3.24498 +0.65000,64.94312,6.79894,38.61004,-3.49243 +0.70000,66.80949,6.61814,36.10108,-3.73957 +0.75000,68.55832,6.42499,33.89831,-3.98644 +0.80000,70.20354,6.21950,31.94888,-4.23310 +0.85000,71.75674,6.00168,30.21148,-4.47958 +0.90000,73.22767,5.77154,28.65330,-4.72590 +0.95000,74.62461,5.52909,27.24796,-4.97210 +1.00000,75.95465,5.27434,25.97403,-5.21818 +1.05000,77.22391,5.00728,24.81390,-5.46417 +1.10000,78.43769,4.72792,23.75297,-5.71007 +1.15000,79.60066,4.43627,22.77904,-5.95590 +1.20000,80.71688,4.13233,21.88184,-6.20166 +1.25000,81.78997,3.81611,21.05263,-6.44737 +1.30000,82.82315,3.48760,20.28398,-6.69302 +1.35000,83.81927,3.14680,19.56947,-6.93863 +1.40000,84.78091,2.79373,18.90359,-7.18420 +1.45000,85.71036,2.42839,18.28154,-7.42973 +1.50000,86.60972,2.05076,17.69912,-7.67522 +1.55000,87.48087,1.66086,17.15266,-7.92069 +1.60000,88.32553,1.25869,16.63894,-8.16612 +1.65000,89.14526,0.84425,16.15509,-8.41153 diff --git a/Verification/Sprinklers_and_Sprays/flat_fire.py b/Verification/Sprinklers_and_Sprays/flat_fire.py new file mode 100644 index 00000000000..b7ced7c78f7 --- /dev/null +++ b/Verification/Sprinklers_and_Sprays/flat_fire.py @@ -0,0 +1,35 @@ + +# This script just computes the exact solution of the flat_fire case. + +import pandas as pd +import numpy as np + +V_0 = 400 +h = 8 +g = 9.8 +C_d = 0.2 +rho_a = 1.2 +rho_d = 1000 +D = 5e-3 + +K = 3 * rho_a * C_d / (4 * rho_d * D) + +dt = 0.05 +tend = 1.65 + +tvec = np.arange(0, tend + dt, dt) + +xexact = np.log(V_0 * K * tvec + 1) / K +yexact = h \ + + (g / (2 * (K * V_0)**2)) * np.log(V_0 * K * tvec + 1) \ + - g * tvec ** 2 / 4 \ + - g * tvec / (2 * K * V_0) +uexact = V_0 / (V_0 * K * tvec + 1) +vexact = g / (2 * K * V_0 * (K * V_0 * tvec + 1)) \ + - g * tvec / 2 \ + - g / (2 * K * V_0) + +df = pd.DataFrame(np.column_stack((tvec, xexact, yexact, uexact, vexact)), columns=['Time', 'x', 'z', 'u', 'w']) + +df.to_csv('flat_fire.csv', float_format='%.5f', index=False) +