Merge pull request #15002 from rmcdermo/master

Manuals: remove MINMOD and MP5 limiter references

Author: rmcdermo

Manuals/FDS_Technical_Reference_Guide/Mass_Chapter.tex

@@ -98,14 +98,12 @@ \subsection{Flux Limiters}
 \hline
 Central Difference                     & 1                              \\
 Godunov                                & 0                              \\
-MINMOD                                 & $\max(0,\min(1,r))$            \\
 Superbee \cite{Roe:1986} (LES default) & $\max(0,\min(2r,1),\min(r,2))$ \\
 CHARM \cite{Zhou:1995} (DNS default)   & $s(3s+1)/(s+1)^2$; $s=1/r$     \\
-MP5 \cite{Suresh:1997}                 & see below
 \end{tabular}
 \end{center}
 \end{table}
-\noindent For the Central Difference, Godunov, MINMOD, and Superbee limiters, the scalar face value is found from
+\noindent For the Central Difference, Godunov, and Superbee limiters, the scalar face value is found from
 \begin{equation}
 \label{eqn_flux_limiter}
 \overline{\phi}^{\rm FL}_{i+1/2} = \left\{ \begin{array}{lcll} \phi_i &+& B(r) \,\frac{1}{2} \,\delta \phi_{\rm loc} & \mbox{if} \quad u_i>0 \vspace{0.2 cm}\\
@@ -117,35 +115,30 @@ \subsection{Flux Limiters}
 \overline{\phi}^{\rm FL}_{i+1/2} = \left\{ \begin{array}{lcll} \phi_i &+& B(r) \,\frac{1}{2} \,\delta \phi_{\rm up} & \mbox{if} \quad u_i>0 \vspace{0.2 cm}\\
 \phi_{i+1} &-& B(r) \,\frac{1}{2} \,\delta \phi_{\rm up} & \mbox{if} \quad u_i<0 \end{array} \right.
 \end{equation}
-The MP5 scheme of Suresh and Huynh \cite{Suresh:1997} is based on the keen observation that three points cannot distinguish between extrema and discontinuities.  The functional form of the limiter is not as simple as the three-point schemes described above, so we refer the reader to the original paper or the FDS source code for details.  But the basic idea behind the method is to use a five-point stencil, three upwind and two downwind, to reconstruct the cell face value, considering both accuracy and monotonicity-preserving constraints.  An additional benefit of the MP5 scheme is that it was designed specifically with strong stability-preserving (SSP) Runge-Kutta time discretizations in mind.  The predictor-corrector scheme used by FDS is similar to the second-order SSP scheme described in \cite{Gottlieb:2001}.
 
 
 \subsubsection{Notes on Implementation}
 
 In practice, we set $r=0$ initially and only compute $r$ if the denominator is not zero.  Note that for $\delta \phi_{loc}=0$, it does not matter which limiter is used: all the limiters yield the same scalar face value.  For CHARM, we set both $r=0$ and $B=0$ initially and only compute $B$ if $r>0$ (this requires data variations to have the same sign). Otherwise, CHARM reduces to Godunov's scheme.
 
-The Central Difference, Godunov, and MINMOD limiters are included for completeness, debugging, and educational purposes.  These schemes have little utility for typical FDS applications.
+The Central Difference and Godunov limiters are included for completeness, debugging, and educational purposes.  These schemes have little utility for typical FDS applications.
 
-\subsubsection{Dealing with Variable Molecular Weights}
+% \subsubsection{Dealing with Variable Molecular Weights}
 
-Maintaining isothermal flow requires
-\begin{equation}
-T = \frac{\overline{W} \bar{p}}{\rho R} = \frac{\bar{p}}{R \rho \sum_\alpha \frac{Z_\alpha}{W_\alpha}} = \frac{\bar{p}}{R \sum_\alpha \frac{\rho Z_\alpha}{W_\alpha}}
-\end{equation}
-to be constant and uniform at all cells and faces. Therefore, with $\bar{p}$ and $R$ constant and uniform, we must maintain
-\begin{equation}
-\label{eq:rho_mw}
-\sum_\alpha \frac{(\rho Z_\alpha)}{W_\alpha} = \frac{\rho}{\overline{W}}
-\end{equation}
+%% leave this here for a moment as a reminder to write up the constant limiter coefficient method we are now using.
 
-The above condition is automatically satisfied in the cases of using Godunov or Central differencing or in the case of binary flow (two species).  However, if we apply a second-order flux limiter, such as Superbee or CHARM, independently to each species in a multi-component (three or more species) flow with variable molecular weights, then this condition is easily violated.
+% Maintaining isothermal flow requires
+% \begin{equation}
+% T = \frac{\overline{W} \bar{p}}{\rho R} = \frac{\bar{p}}{R \rho \sum_\alpha \frac{Z_\alpha}{W_\alpha}} = \frac{\bar{p}}{R \sum_\alpha \frac{\rho Z_\alpha}{W_\alpha}}
+% \end{equation}
+% to be constant and uniform at all cells and faces. Therefore, with $\bar{p}$ and $R$ constant and uniform, we must maintain
+% \begin{equation}
+% \label{eq:rho_mw}
+% \sum_\alpha \frac{(\rho Z_\alpha)}{W_\alpha} = \frac{\rho}{\overline{W}}
+% \end{equation}
+
+% The above condition is automatically satisfied in the cases of using Godunov or Central differencing or in the case of binary flow (two species).  However, if we apply a second-order flux limiter, such as Superbee or CHARM, independently to each species in a multi-component (three or more species) flow with variable molecular weights, then this condition is easily violated.
 
-To handle this situation, in LES mode, FDS will apply a correction to the most abundant species locally.  We first compute the flux-limited face values of the mass density over the mixture-average molecular weight.  Then we compute flux-limited face values of the species densities.  Finally, the error in Eq.~(\ref{eq:rho_mw}) is absorbed into the most abundant species locally,
-\begin{equation}
-\label{eq:rhoZ_cor}
-\overline{\rho Z_\alpha}^{\rm COR} = W_\alpha \left( \overline{\left\{\frac{\rho}{\overline{W}}\right\}}^{\rm FL} - \sum_{\beta \ne \alpha} \frac{\overline{\{\rho Z_\beta\}}^{\rm FL}}{W_\beta} \right)
-\end{equation}
-where $\alpha$ is the most abundant species on the face.
 
 \subsection{Time Splitting for Mass Source Terms}
 \label{sec_time_splitting}

Manuals/FDS_User_Guide/FDS_User_Guide.tex

@@ -9123,7 +9123,7 @@ \subsection{Heat Transfer Constraint}
 \section{Flux Limiters}
 \label{info:flux_limiters}
 
-FDS employs \emph{total variation diminishing} (TVD) schemes for scalar transport.  The default for VLES (FDS default \ct{SIMULATION_MODE}) is Superbee \cite{Roe:1986}, so chosen because this scheme does the best job preserving the scalar variance in highly turbulent flows with coarse grid resolution.  The default scheme for DNS and LES is CHARM \cite{Zhou:1995} because the gradient steepening used in Superbee forces a stair step pattern at high resolution, while CHARM is convergent.   A few other schemes (including Godunov and central differencing) are included for completeness; more details can be found in the Tech Guide \cite{FDS_Tech_Guide}.  Table \ref{tab:flux_limiters} below shows the character strings which may be used to invoke the various limiter schemes.
+FDS employs \emph{total variation diminishing} (TVD) schemes for scalar transport.  The default for VLES (FDS default \ct{SIMULATION_MODE}) is Superbee \cite{Roe:1986}, so chosen because this scheme does the best job preserving the scalar variance in highly turbulent flows with coarse grid resolution.  The default scheme for DNS and LES is CHARM \cite{Zhou:1995} because the gradient steepening used in Superbee forces a stair step pattern at high resolution, while CHARM is convergent.   Godunov and central differencing are included for completeness; more details can be found in the Tech Guide \cite{FDS_Tech_Guide}.  Table \ref{tab:flux_limiters} below shows the character strings which may be used to invoke the various limiter schemes.
 
 \begin{lstlisting}
 &MISC FLUX_LIMITER='GODUNOV' / ! invoke Godunov (first-order upwind scheme)
@@ -9140,9 +9140,7 @@ \section{Flux Limiters}
 Central differencing                & \ct{'CENTRAL'}  \\
 Godunov                             & \ct{'GODUNOV'}  \\
 Superbee (VLES, SVLES default)      & \ct{'SUPERBEE'}  \\
-MINMOD                              & \ct{'MINMOD'}  \\
 CHARM (DNS, LES default)            & \ct{'CHARM'}  \\
-MP5                                 & \ct{'MP5'}  \\ \hline
 \end{tabular}
 \end{table}
 

Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex

@@ -893,16 +893,15 @@ \subsection{Solid Body Rotation Velocity Field (\texorpdfstring{\ct{soborot}}{so
    \label{fig_soborot_square_wave}
 \end{figure}
 
-For the cosine wave initial condition the derivatives of the scalar field are continuous.  Therefore, we see ${\cal O}(\delta x^2)$ convergence of the CHARM and MP5 schemes, as shown in Fig.~\ref{fig_soborot_cos_wave}.  Superbee shows smaller error at coarse resolution, but the gradient steepening degenerates its accuracy at higher resolutions---hence CHARM is selected for LES and DNS.
+For the cosine wave initial condition the derivatives of the scalar field are continuous.  Therefore, we see ${\cal O}(\delta x^2)$ convergence of the CHARM scheme, as shown in Fig.~\ref{fig_soborot_cos_wave}.  Superbee shows smaller error at coarse resolution, but the gradient steepening degenerates its accuracy at higher resolutions---hence CHARM is selected for LES and DNS.
 
 \begin{figure}[ht]
    \begin{tabular*}{\textwidth}{l@{\extracolsep{\fill}}r}
       \includegraphics[height=2.2in]{SCRIPT_FIGURES/soborot_superbee_cos_wave} &
       \includegraphics[height=2.2in]{SCRIPT_FIGURES/soborot_charm_cos_wave}   \\
-      \includegraphics[height=2.2in]{SCRIPT_FIGURES/soborot_mp5_cos_wave} &
-      \includegraphics[height=2.2in]{SCRIPT_FIGURES/soborot_cos_wave_error}
+      \includegraphics[height=2.2in]{SCRIPT_FIGURES/soborot_cos_wave_error} &
    \end{tabular*}
-   \caption[Solid body rotation cosine wave convergence]{Solid body rotation cosine wave solution convergence.  Scalar fields along upper-left diagonal for Superbee (upper-left), CHARM (upper-right), and MP5 (lower-left).  (Lower-right) L2 Error.}
+   \caption[Solid body rotation cosine wave convergence]{Solid body rotation cosine wave solution convergence.  Scalar fields along upper-left diagonal for Superbee (upper-left) and CHARM (upper-right).  (Lower-left) L2 Error.}
    \label{fig_soborot_cos_wave}
 \end{figure}
 

Utilities/Matlab/scripts/compression_wave.m

@@ -61,23 +61,6 @@
     skip_case = 1;
 end
 
-if ~exist([data_dir,'compression_wave_FL5_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_FL5_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_FL5_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_FL5_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
@@ -113,16 +96,6 @@
 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);
 
-% MP5 limiter, FL=5
-M_FL5_16 = csvread([data_dir,'compression_wave_FL5_16_devc.csv'],2,0);
-M_FL5_32 = csvread([data_dir,'compression_wave_FL5_32_devc.csv'],2,0);
-M_FL5_64 = csvread([data_dir,'compression_wave_FL5_64_devc.csv'],2,0);
-M_FL5_128 = csvread([data_dir,'compression_wave_FL5_128_devc.csv'],2,0);
-t_FL5_16 = M_FL5_16(:,1); rho_fds_FL5_16 = M_FL5_16(:,3);
-t_FL5_32 = M_FL5_32(:,1); rho_fds_FL5_32 = M_FL5_32(:,3);
-t_FL5_64 = M_FL5_64(:,1); rho_fds_FL5_64 = M_FL5_64(:,3);
-t_FL5_128 = M_FL5_128(:,1); rho_fds_FL5_128 = M_FL5_128(:,3);
-
 % analytical solution
 
 a = 2;
@@ -147,11 +120,6 @@
 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);
 
-rho_FL5_16 = compression_wave_soln(rho_fds_FL5_16(1),x-L/32,y-L/32,a,c,t_FL5_16);      error_FL5_16 = norm(rho_fds_FL5_16-rho_FL5_16)/length(t_FL5_16);
-rho_FL5_32 = compression_wave_soln(rho_fds_FL5_32(1),x-L/64,y-L/64,a,c,t_FL5_32);      error_FL5_32 = norm(rho_fds_FL5_32-rho_FL5_32)/length(t_FL5_32);
-rho_FL5_64 = compression_wave_soln(rho_fds_FL5_64(1),x-L/128,y-L/128,a,c,t_FL5_64);    error_FL5_64 = norm(rho_fds_FL5_64-rho_FL5_64)/length(t_FL5_64);
-rho_FL5_128 = compression_wave_soln(rho_fds_FL5_128(1),x-L/256,y-L/256,a,c,t_FL5_128); error_FL5_128 = norm(rho_fds_FL5_128-rho_FL5_128)/length(t_FL5_128);
-
 figure
 plot_style
 set(gca,'Units',Plot_Units)
@@ -197,20 +165,18 @@
 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];
-e_FL5 = [error_FL5_16 error_FL5_32 error_FL5_64 error_FL5_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,e_FL5,'c>-','LineWidth',Line_Width); hold on
-H(5)=loglog(h,.1*h,'k--','LineWidth',Line_Width);
-H(6)=loglog(h,.1*h.^2,'k-','LineWidth',Line_Width);
+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:6),'Central','Superbee','CHARM','MP5','{\it O}({\it\deltax})','{\it O}({\it\deltax^2})','Location','NorthWest');
+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');

Utilities/Matlab/scripts/soborot_mass_transport.m

@@ -23,10 +23,6 @@
         'soborot_godunov_square_wave_16', ...
         'soborot_godunov_square_wave_32', ...
         'soborot_godunov_square_wave_64', ...
-        'soborot_mp5_cos_wave_128', ...
-        'soborot_mp5_cos_wave_16', ...
-        'soborot_mp5_cos_wave_32', ...
-        'soborot_mp5_cos_wave_64', ...
         'soborot_superbee_cos_wave_128', ...
         'soborot_superbee_cos_wave_16', ...
         'soborot_superbee_cos_wave_32', ...
@@ -423,62 +419,6 @@
 set(gcf,'Position',[0 0 Paper_Width Paper_Height]);
 print(gcf,'-dpdf',[plot_dir,'soborot_superbee_cos_wave']);
 
-% FLUX_LIMITER='MP5'
-
-M = importdata([data_dir,'soborot_mp5_cos_wave_16_devc.csv'],',',2);
-col_start = find(strcmp(M.colheaders,'"Y_TRACER-1"'));
-col_end   = find(strcmp(M.colheaders,'"Y_TRACER-16"'));
-Y_mp5_16 = M.data(end,col_start:col_end);
-
-M = importdata([data_dir,'soborot_mp5_cos_wave_32_devc.csv'],',',2);
-col_start = find(strcmp(M.colheaders,'"Y_TRACER-1"'));
-col_end   = find(strcmp(M.colheaders,'"Y_TRACER-32"'));
-Y_mp5_32 = M.data(end,col_start:col_end);
-
-M = importdata([data_dir,'soborot_mp5_cos_wave_64_devc.csv'],',',2);
-col_start = find(strcmp(M.colheaders,'"Y_TRACER-1"'));
-col_end   = find(strcmp(M.colheaders,'"Y_TRACER-64"'));
-Y_mp5_64 = M.data(end,col_start:col_end);
-
-M = importdata([data_dir,'soborot_mp5_cos_wave_128_devc.csv'],',',2);
-col_start = find(strcmp(M.colheaders,'"Y_TRACER-1"'));
-col_end   = find(strcmp(M.colheaders,'"Y_TRACER-128"'));
-Y_mp5_128 = M.data(end,col_start:col_end);
-
-figure
-set(gca,'Units',Plot_Units)
-set(gca,'Position',[Plot_X Plot_Y Plot_Width Plot_Height])
-Y_exact = zeros(1,length(r_exact));
-i_range = find(r_exact>0.25 & r_exact<0.75);
-Y_exact(i_range) = 0.5*(1. + cos(4.*pi*r_exact(i_range))); % see PERIODIC_TEST==13 in wall.f90
-H(1)=plot(r_exact,Y_exact,'k-'); hold on
-H(2)=plot(r_16,Y_mp5_16,'bo-');
-H(3)=plot(r_32,Y_mp5_32,'m*-');
-H(4)=plot(r_64,Y_mp5_64,'r^-');
-H(5)=plot(r_128,Y_mp5_128,'gsq-');
-
-xlabel('Radial Position (m)','Interpreter',Font_Interpreter,'FontSize',Label_Font_Size)
-ylabel('Scalar Mass Fraction','Interpreter',Font_Interpreter,'FontSize',Label_Font_Size)
-axis([0 1 0 1.2])
-text(.05,1.1,'MP5','FontName',Font_Name,'FontSize',Label_Font_Size)
-
-set(gca,'FontName',Font_Name)
-set(gca,'FontSize',Label_Font_Size)
-
-lh=legend(H,'Exact','{\itn}=16','{\itn}=32','{\itn}=64','{\itn}=128');
-set(lh,'FontName',Font_Name,'FontSize',Key_Font_Size)
-legend boxoff
-
-Git_Filename = [data_dir,'soborot_mp5_cos_wave_16_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',[plot_dir,'soborot_mp5_cos_wave']);
-
 % plot error
 
 Y_exact_16 = interp1(r_exact,Y_exact,r_16);
@@ -496,11 +436,6 @@
 e_superbee_64 = norm(Y_superbee_64-Y_exact_64,1)/length(r_64);
 e_superbee_128 = norm(Y_superbee_128-Y_exact_128,1)/length(r_128);
 
-e_mp5_16 = norm(Y_mp5_16-Y_exact_16,1)/length(r_16);
-e_mp5_32 = norm(Y_mp5_32-Y_exact_32,1)/length(r_32);
-e_mp5_64 = norm(Y_mp5_64-Y_exact_64,1)/length(r_64);
-e_mp5_128 = norm(Y_mp5_128-Y_exact_128,1)/length(r_128);
-
 dx = L./[16 32 64 128];
 
 figure
@@ -510,7 +445,6 @@
 H(2)=loglog(dx,dx.^2,'k--');
 H(3)=loglog(dx,[e_charm_16 e_charm_32 e_charm_64 e_charm_128],'ko-');
 H(4)=loglog(dx,[e_superbee_16 e_superbee_32 e_superbee_64 e_superbee_128],'k*-');
-H(5)=loglog(dx,[e_mp5_16 e_mp5_32 e_mp5_64 e_mp5_128],'ksq-');
 
 % trap cosine wave error
 if e_charm_128>2.9e-04
@@ -519,9 +453,6 @@
 if e_superbee_128>6.4e-04
     display(['Error: soborot_charm_square_wave_128 out of tolerance'])
 end
-if e_mp5_128>7.9e-05
-    display(['Error: soborot_superbee_square_wave_128 out of tolerance'])
-end
 
 axis([min(dx) .1 min(dx.^2) max(dx)])
 xlabel('Grid Spacing (m)','Interpreter',Font_Interpreter,'FontSize',Label_Font_Size)
@@ -531,7 +462,7 @@
 set(gca,'FontName',Font_Name)
 set(gca,'FontSize',Label_Font_Size)
 
-lh=legend(H,'{\itO}(\delta{\itx})','{\itO}(\delta{\itx}^2)','CHARM','Superbee','MP5');
+lh=legend(H,'{\itO}(\delta{\itx})','{\itO}(\delta{\itx}^2)','CHARM','Superbee');
 set(lh,'Location','Southeast')
 set(lh,'FontName',Font_Name,'FontSize',Key_Font_Size)
 legend boxoff