source: anuga_core/documentation/user_manual/old_pyvolution_documentation/limiting.tex @ 4265

\documentclass[12pt]{article}
\usepackage{graphicx}


% THIS HAS BEEN MOVED INTO THE USER MANUAL (16 Feb 2007)

\newcommand{\hmin}{h_{\mbox{\tiny min}}}

\begin{document}

Ole's thoughts on flux and slope limiting in ANUGA - 7 February 2007.

This is a rough draft, please let me know if this makes sense.
O

\section{Flux limiting}

The shallow water equations are solved numerically using a 
finite volume method on unstructured triangular grid.
The upwind central scheme due to Kurganov and Petrova is used as an 
approximate Riemann solver for the computation of inviscid flux functions.
This makes it possible to handle discontinuous solutions. 

To alleviate the problems associated with numerical instabilities due to 
small water depths near a wet/dry boundary we employ a new flux limiter that
ensures that unphysical fluxes are never encounted.


Let $u$ and $v$ be the velocity components in the $x$ and $y$ direction, 
$w$ the absolute water level (stage) and
$z$ the bed elevation. The latter are assumed to be relative to the 
same height datum. 
The conserved quantities tracked by ANUGA are momentum in the 
$x$-direction ($\mu = uh$), momentum in the $y$-direction ($\nu = vh$) 
and depth ($h = w-z$).

The flux calculation requires access to the velocity vector $(u, v)$ 
where each component is obtained as $u = \mu/h$ and $v = \nu/h$ respectively.
In the presence of very small water depths, these calculations become 
numerically unreliable and will typically cause unphysical speeds.

We have employed a flux limiter which replaces the calculations above with 
the limited approximations.
\begin{equation}
  \hat{u} = \frac{\mu}{h + h_0/h}, \bigskip \hat{v} = \frac{\nu}{h + h_0/h}, 
\end{equation}
where $h_0$ is a regularisation parameter that controls the minimal 
magnitude of the denominator. Taking the limits we have for $\hat{u}$
\[
  \lim_{h \rightarrow 0} \hat{u} = 
  \lim_{h \rightarrow 0} \frac{\mu}{h + h_0/h} = 0
\]
and 
\[
  \lim_{h \rightarrow \infty} \hat{u} = 
  \lim_{h \rightarrow \infty} \frac{\mu}{h + h_0/h} = \frac{\mu}{h} = u
\]
with similar results for $\hat{v}$.

The maximal value of $\hat{u}$ is attained when $h+h_0/h$ is minimal or (by differentiating the denominator)
\[
  1 - h_0/h^2 = 0
\]
or
\[
  h_0 = h^2
\]


ANUGA has a global parameter $H_0$ that controls the minimal depth which 
is considered in the various equations. This parameter is typically set to
$10^{-3}$. Setting
\[
  h_0 = H_0^2
\]
provides a reasonable balance between accurracy and stability. In fact,
setting $h=H_0$ will scale the predicted speed by a factor of $0.5$:
\[
  \left[ \frac{\mu}{h + h_0/h} \right]_{h = H_0} = \frac{\mu}{2 H_0}
\]
In general, for multiples of the minimal depth $N H_0$ one obtains 
\[
  \left[ \frac{\mu}{h + h_0/h} \right]_{h = N H_0} = 
  \frac{\mu}{H_0 (1 + 1/N^2)}
\]
which converges quadratically to the true value with the multiple N. 


%The developed numerical model has been applied to several test cases as well as to real flows. Numerical tests prove the robustness and accuracy of the model.





\section{Slope limiting}
A multidimensional slope-limiting technique is employed to achieve second-order spatial accuracy and to prevent spurious oscillations. This is using the MinMod limiter and is documented elsewhere.

However close to the bed, the limiter must ensure that no negative depths occur. On the other hand, in deep water, the bed topography should be ignored for the purpose of the limiter.


Let $w, z, h$  be the stage, bed elevation and depth at the centroid and
let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$.
Define the minimal depth across all vertices as $\hmin$ as
\[
  \hmin = \min_i h_i
\]

Let $\tilde{w_i}$ be the stage obtained from a gradient limiter 
limiting on stage only. The corresponding depth is then defined as
\[
  \tilde{h_i} = \tilde{w_i} - z_i
\]
We would use this limiter in deep water which we will define (somewhat boldly) 
as
\[
  \hmin \ge \epsilon
\]


Similarly, let $\bar{w_i}$ be the stage obtained from a gradient 
limiter limiting on depth respecting the bed slope.
The corresponding depth is defined as
\[
  \bar{h_i} = \bar{w_i} - z_i
\]


We introduce the concept of a balanced stage $w_i$ which is obtained as
the linear combination

\[
  w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i}
\]
or
\[
  w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i}
\]
where $\alpha \in [0, 1]$.

Since $\tilde{w_i}$ is obtained in 'deep' water where the bedslope 
is ignored we have immediately that 
\[
  \alpha = 1 \mbox{ for } \hmin \ge \epsilon %or dz=0
\]
%where the maximal bed elevation range $dz$ is defined as
%\[
%  dz = \max_i |z_i - z|
%\]

If $\hmin < \epsilon$ we want to use the 'shallow' limiter just enough that
no negative depths occur. Formally, we will require that
\[
  \alpha \tilde{h_i} + (1-\alpha) \bar{h_i} > \epsilon, \forall i
\]
or 
\begin{equation} 
  \alpha(\tilde{h_i} - \bar{h_i}) > \epsilon - \bar{h_i}, \forall i
  \label{eq:limiter bound}
\end{equation} 

There are two cases:
\begin{enumerate} 
  \item $\bar{h_i} \le \tilde{h_i}$: The deep water (limited using stage) 
  vertex is at least as far away from the bed than the shallow water 
  (limited using depth). In this case we won't need any contribution from
  $\bar{h_i}$ and can accept any $alpha$. 
  
  E.g.\ $\alpha=1$ reduces Equation \ref{eq:limiter bound} to
  \[
    \tilde{h_i} > \epsilon  
  \]
  whereas $\alpha=0$ yields
  \[
    \bar{h_i} > \epsilon        
  \]
  all well and good.
  \item $\bar{h_i} > \tilde{h_i}$: In this case the the deep water vertex is 
  closer to the bed than the shallow water vertex or even below the bed. 
  In this case we need to find an $alpha$ that will ensure a positive depth.   
  Rearranging Equation \ref{eq:limiter bound} and solving for $\alpha$ one 
  obtains the bound
  \[
    \alpha < \frac{\epsilon - \bar{h_i}}{\tilde{h_i} - \bar{h_i}}, \forall i
  \] 
\end{enumerate} 

Ensuring Equation \ref{eq:limiter bound} holds true for all vertices one 
arrives at the definition
\[
  \alpha = \min_{i} \frac{\bar{h_i} - \epsilon}{\bar{h_i} - \tilde{h_i}}
\]
which will guarantee that no vertex 'cuts' through the bed. Finally, should
$\bar{h_i} < \epsilon$ and therefore $\alpha < 0$, we suggest setting 
$alpha=0$ and similarly capping $\alpha$ at 1 just in case.

%Furthermore,
%dropping the $\epsilon$ ensures that alpha is always positive and also  
%provides a numerical safety {??)
  
  









\section{Slope limiting (old version)}

Let $w, z, h$  be the stage, bed elevation and depth at the centroid and
let $w_i, z_i, h_i$ be the stage, bed elevation and depth at vertex $i$.

Define the maximal bed elevation range $dz$ as

\[
  dz = \max_i |z_i - z|
\]

and the minimal depth $\hmin$ as

\[
  \hmin = \min_i h_i
\]


\[
  \alpha = \left \{
  \begin{array}{ll} 
    \max (\min ( 2 \hmin / dz )) & dz > 0 \\
    1 & dz \leq 0     
  \end{array}
  \right .  
\]


Let $\tilde{w_i}$ be the stage obtained from a gradient limiter limiting on stage. The corresponding depth is the defined as

\[
  \tilde{h_i} = \tilde{w_i} - z_i
\]

Let $\bar{h_i}$ be the depth obtained from a gradient limiter limiting on depth.
The corresponding stage is the defined as

\[
  \bar{w_i} = z_i + \bar{h_i}
\]


The balanced stage $w_i$ is then obtained by the linear combination

\[
  w_i = \alpha \tilde{w_i} + (1-\alpha) \bar{w_i}
\]

or

\[
  w_i = z_i + \alpha \tilde{h_i} + (1-\alpha) \bar{h_i}
\]

\end{document}