Skip to content

Manuals: remove MINMOD and MP5 limiter references #15002

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
rmcdermo merged 3 commits into from
Aug 21, 2025
Merged

Manuals: remove MINMOD and MP5 limiter references #15002

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5e1b051
Matlab: remove MP5 scheme verification scripts
rmcdermo Aug 21, 2025
140d661
Manuals: remove MINMOD and MP5
rmcdermo Aug 21, 2025
9b7ca57
FDS Tech Guide: remove refs to MP5 and MINMOD and MW correction
rmcdermo Aug 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 15 additions & 22 deletions Manuals/FDS_Technical_Reference_Guide/Mass_Chapter.tex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}\\
Expand All @@ -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}
Expand Down
4 changes: 1 addition & 3 deletions Manuals/FDS_User_Guide/FDS_User_Guide.tex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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}

Expand Down
7 changes: 3 additions & 4 deletions Manuals/FDS_Verification_Guide/FDS_Verification_Guide.tex
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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}

Expand Down
40 changes: 3 additions & 37 deletions Utilities/Matlab/scripts/compression_wave.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand All @@ -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)
Expand Down Expand Up @@ -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');
Expand Down
71 changes: 1 addition & 70 deletions Utilities/Matlab/scripts/soborot_mass_transport.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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', ...
Expand Down Expand Up @@ -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);
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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)
Expand All @@ -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
Expand Down