Changeset 4265

Timestamp:

Feb 16, 2007, 5:26:04 PM (17 years ago)

Author:

ole

Message:

Move mathematical destcription from linuxconf paper and limiting.tex into user manual.
Have not tried to compile it though.

Location:

anuga_core/documentation/user_manual

Files:

• anuga_user_manual.tex (modified) (1 diff)
• definitions.tex (modified) (1 diff)
• old_pyvolution_documentation/limiting.tex (modified) (1 diff)

anuga_core/documentation/user_manual/anuga_user_manual.tex

@@ -2660 +2660 @@
   Elevation Model) and converts it to PTS format.
   \end{funcdesc}
+
+
+  
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\chapter{\anuga mathematical background}
+\label{cd:mathematical background}
+
+\section{Introduction}
+
+This chapter outlines the mathematics underpinning \anuga.
+  
+  
+
+\section{Model}
+\label{sec:model}
+
+The shallow water wave equations are a system of differential
+conservation equations which describe the flow of a thin layer of
+fluid over terrain. The form of the equations are:
+\[
+\frac{\partial \UU}{\partial t}+\frac{\partial \EE}{\partial
+x}+\frac{\partial \GG}{\partial y}=\SSS
+\]
+where $\UU=\left[ {{\begin{array}{*{20}c}
+ h & {uh} & {vh} \\
+\end{array} }} \right]^T$ is the vector of conserved quantities; water depth
+$h$, $x$-momentum $uh$ and $y$-momentum $vh$. Other quantities
+entering the system are bed elevation $z$ and stage (absolute water
+level) $w$, where the relation $w = z + h$ holds true at all times.
+The fluxes in the $x$ and $y$ directions, $\EE$ and $\GG$ are given
+by
+\[
+\EE=\left[ {{\begin{array}{*{20}c}
+ {uh} \hfill \\
+ {u^2h+gh^2/2} \hfill \\
+ {uvh} \hfill \\
+\end{array} }} \right]\mbox{ and }\GG=\left[ {{\begin{array}{*{20}c}
+ {vh} \hfill \\
+ {vuh} \hfill \\
+ {v^2h+gh^2/2} \hfill \\
+\end{array} }} \right]
+\]
+and the source term (which includes gravity and friction) is given
+by
+\[
+\SSS=\left[ {{\begin{array}{*{20}c}
+ 0 \hfill \\
+ -{gh(z_{x} + S_{fx} )} \hfill \\
+ -{gh(z_{y} + S_{fy} )} \hfill \\
+\end{array} }} \right]
+\]
+where $S_f$ is the bed friction. The friction term is modelled using
+Manning's resistance law
+\[
+S_{fx} =\frac{u\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}\mbox{ and }S_{fy}
+=\frac{v\eta ^2\sqrt {u^2+v^2} }{h^{4/3}}
+\]
+in which $\eta$ is the Manning resistance coefficient.
+
+As demonstrated in our papers, \cite{modsim2005,Rob99l} these
+equations provide an excellent model of flows associated with
+inundation such as dam breaks and tsunamis.
+
+\section{Finite Volume Method}
+\label{sec:fvm}
+
+We use a finite-volume method for solving the shallow water wave
+equations \cite{Rob99l}. The study area is represented by a mesh of
+triangular cells as in Figure~\ref{fig:mesh} in which the conserved
+quantities of  water depth $h$, and horizontal momentum $(uh, vh)$,
+in each volume are to be determined. The size of the triangles may
+be varied within the mesh to allow greater resolution in regions of
+particular interest.
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=5.0cm,keepaspectratio=true]{step-five}
+\caption{Triangular mesh used in our finite volume method. Conserved
+quantities $h$, $uh$ and $vh$ are associated with the centroid of
+each triangular cell.} \label{fig:mesh}
+\end{center}
+\end{figure}
+
+The equations constituting the finite-volume method are obtained by
+integrating the differential conservation equations over each
+triangular cell of the mesh. Introducing some notation we use $i$ to
+refer to the $i$th triangular cell $T_i$, and ${\cal N}(i)$ to the
+set of indices referring to the cells neighbouring the $i$th cell.
+Then $A_i$ is the area of the $i$th triangular cell and $l_{ij}$ is
+the length of the edge between the $i$th and $j$th cells.
+
+By applying the divergence theorem we obtain for each volume an
+equation which describes the rate of change of the average of the
+conserved quantities within each cell, in terms of the fluxes across
+the edges of the cells and the effect of the source terms. In
+particular, rate equations associated with each cell have the form
+$$
+ \frac{d\UU_i }{dt}+ \frac1{A_i}\sum\limits_{j\in{\cal N}(i)} \HH_{ij} l_{ij} = \SSS_i
+$$
+where
+\begin{itemize}
+\item $\UU_i$ the vector of conserved quantities averaged over the $i$th cell,
+\item $\SSS_i$ is the source term associated with the $i$th cell,
+and
+\item $\HH_{ij}$ is the outward normal flux of
+material across the \textit{ij}th edge.
+\end{itemize}
+
+
+%\item $l_{ij}$ is the length of the edge between the $i$th and $j$th
+%cells
+%\item $m_{ij}$ is the midpoint of
+%the \textit{ij}th edge,
+%\item
+%$\mathbf{n}_{ij} = (n_{ij,1} , n_{ij,2})$is the outward pointing
+%normal along the \textit{ij}th edge, and The
+
+The flux $\HH_{ij}$ is evaluated using a numerical flux function
+$\HH(\cdot, \cdot ; \ \cdot)$ which is consistent with the shallow
+water flux in the sense that for all conservation vectors $\UU$ and normal vectors $\nn$
+$$
+H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 .
+$$
+
+Then
+$$
+\HH_{ij}  = \HH(\UU_i(m_{ij}),
+\UU_j(m_{ij}); \mathbf{n}_{ij})
+$$
+where $m_{ij}$ is the midpoint of the \textit{ij}th edge and
+$\mathbf{n}_{ij}$ is the outward pointing normal, with respect to the $i$th cell, on the
+\textit{ij}th edge. The function $\UU_i(x)$ for $x \in
+T_i$ is obtained from the vector $\UU_k$ of conserved average values for the $i$th and
+neighbouring  cells.
+
+We use a second order reconstruction to produce a piece-wise linear
+function construction of the conserved quantities for  all $x \in
+T_i$ for each cell (see Figure~\ref{fig:mesh:reconstruct}. This
+function is allowed to be discontinuous across the edges of the
+cells, but the slope of this function is limited to avoid
+artificially introduced oscillations.
+
+Godunov's method (see \cite{Toro-92}) involves calculating the
+numerical flux function $\HH(\cdot, \cdot ; \ \cdot)$ by exactly
+solving the corresponding one dimensional Riemann problem normal to
+the edge. We use the central-upwind scheme of \cite{KurNP2001} to
+calculate an approximation of the flux across each edge.
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=5.0cm,keepaspectratio=true]{step-reconstruct}
+\caption{From the values of the conserved quantities at the centroid
+of the cell and its neighbouring cells, a discontinuous piecewise
+linear reconstruction of the conserved quantities is obtained.}
+\label{fig:mesh:reconstruct}
+\end{center}
+\end{figure}
+
+In the computations presented in this paper we use an explicit Euler
+time stepping method with variable timestepping adapted to the
+observed CFL condition.
+
+
+\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 {??)
+  
+  
 
 

anuga_core/documentation/user_manual/definitions.tex

@@ -10 +10 @@
 %\newcommand{\anugav}[1]{\textbf{ANUGA}$^\copyright$ V#1\ }
 
+
+\newcommand{\Python}{\textsc{Python}}
+\newcommand{\VPython}{\textsc{VPython}}
+\newcommand{\pypar}{\textsc{mpi}}
+\newcommand{\Metis}{\textsc{Metis}}
+\newcommand{\mpi}{\textsc{mpi}}
+
+\newcommand{\UU}{\mathbf{U}}
+\newcommand{\VV}{\mathbf{V}}
+\newcommand{\EE}{\mathbf{E}}
+\newcommand{\GG}{\mathbf{G}}
+\newcommand{\FF}{\mathbf{F}}
+\newcommand{\HH}{\mathbf{H}}
+\newcommand{\SSS}{\mathbf{S}}
+\newcommand{\nn}{\mathbf{n}}
+
+%\newcommand{\code}[1]{\texttt{#1}}
+
+
+
+
+
+
 %Used with \code.
 %

anuga_core/documentation/user_manual/old_pyvolution_documentation/limiting.tex

@@ +1 @@
+
 \documentclass[12pt]{article}
 \usepackage{graphicx}
 
+
+% THIS HAS BEEN MOVED INTO THE USER MANUAL (16 Feb 2007)
 
 \newcommand{\hmin}{h_{\mbox{\tiny min}}}