source: anuga_work/publications/linuxconf_2006/paper_ole_nielsen.tex @ 4115

% Use the standard \LaTeXe\ article style in 12pt Computer Modern font
% on A4 paper by
\documentclass[12pt,a4paper]{article}
% Do \emph{not} change the width nor the height of the text from the
% defaults set by this document class.
%
% The alternative which is closer to what we actually use is
% \documentclass[11pt,a5paper]{article}
% \usepackage[a5paper]{geometry}
% Because it is a great size for on screen reading
% and prints nicely on a4paper either 2up or booklet.

% The preamble is to contain your own \LaTeX\ commands and to say
% what packages to use.  Three useful packages are the following:
\usepackage{graphicx} % avoid epsfig or earlier such packages
\usepackage{url}      % for URLs and DOIs
\newcommand{\doi}[1]{\url{http://dx.doi.org/#1}}
\usepackage{amsmath}  % many want amsmath extensions
\usepackage{amsfonts}
\usepackage{underscore}
% Avoid loading unused packages (as done by some \LaTeX\ editors).

% Create title and authors using \verb|\maketitle|.  Separate authors by
% \verb|\and| and put addresses in \verb|\thanks| command with
% \verb|\url| command \verb|\protect|ed.
\title{Open Source Software for Computational Modelling of Fluid Flow}

\author{
O.~M.~Nielsen\thanks{Risk Assessment Methods Project, Geospatial and
Earth Monitoring Division, Geoscience Australia, Symonston,
\textsc{Australia}. \protect\url{mailto:Ole.Nielsen@ga.gov.au}} 
\and
S.~G.~Roberts\thanks{Dept. of Maths, Australian National University,
Canberra, \textsc{Australia}. \protect\url{mailto:stephen.roberts}}}

\date{30 December 2006}

\newcommand{\AnuGA}{\textsc{Anuga}}
\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}}


%\graphicspath{{../Figures/}}

\begin{document}

% Use default \verb|\maketitle|.
\maketitle

% Use the \verb|abstract| environment.
\begin{abstract}
Geoscience Australia and the Australian National University are
developing a hydrodynamic inundation modelling tool called \AnuGA{}
to help simulate the impact of natural inundation hazards such as
riverine flooding, storm surges and tsunami. The core of \AnuGA{} is
a \Python{} implementation of a finite-volume method for solving the
conservative form of the Shallow Water Wave equation.  In this paper
we describe the parallelisation of the code using a domain
decomposition strategy. We describe the use of the the \Metis{}
graph partitioning library to decompose our finite volume meshes.
The parallel efficiency of our code is tested using a number of mesh
partitions, and we verify that the \Metis{} graph partition is
particularly efficient.
\end{abstract}

% By default we include a table of contents in each paper.
\tableofcontents

% Use \verb|\section|, \verb|\subsection|, \verb|\subsubsection| and
% possibly \verb|\paragraph| to structure your document.  Make sure
% you \verb|\label| them for cross-referencing with \verb|\ref|\,.
\section{Introduction}
\label{sec:intro}

%Floods are the single greatest cause of death due to natural hazards
%in Australia, causing almost 40{\%} of the fatalities recorded
%between 1788 and 2003~\cite{Blong-2005}. Analysis of buildings
%damaged between 1900 and 2003 suggests that 93.6{\%} of damage is
%the result of meteorological hazards, of which almost 25{\%} is
%directly attributable to flooding~\cite{Blong-2005}.

Flooding of coastal communities may result from surges of near-shore
waters caused by severe storms. The extent of inundation is
critically linked to tidal conditions, bathymetry and topography; as
recently exemplified in the United States by Hurricane Katrina.
While the scale of the impact from such events is not common, the
preferential development of Australian coastal corridors means that
storm-surge inundation of even a few hundred metres beyond the
shoreline has increased potential to cause significant disruption
and loss.

Coastal communities also face the small but real risk of tsunami.
Fortunately, catastrophic tsunami of the scale of the 26 December
2004 event are exceedingly rare. However, smaller-scale tsunami are
more common and regularly threaten coastal communities around the
world. Earthquakes which occur in the Java Trench near Indonesia
(e.g.~\cite{TsuMIS1995}) and along the Puysegur Ridge to the south
of New Zealand (e.g.~\cite{LebKC1998}) have potential to generate
tsunami that may threaten Australia's northwestern and southeastern
coastlines.

Hydrodynamic modelling allows flooding, storm-surge and tsunami
hazards to be better understood, their impacts to be anticipated
and, with appropriate planning, their effects to be mitigated.
Geoscience Australia in collaboration with the Mathematical Sciences
Institute, Australian National University, is developing a software
application called \AnuGA{} to model the hydrodynamics of floods,
storm surges and tsunami. These hazards are modelled using the
conservative shallow water equations which are described in
section~\ref{sec:model}. In \AnuGA{} these equations are solved
using a finite volume method as described in section~\ref{sec:fvm}.
A more complete discussion of the method can be found in
\cite{modsim2005} where the model and solution technique is
validated on a standard tsunami benchmark data set.

In this paper we will provide a description of the parallelisation
of the code using a domain decomposition strategy and provide the
preliminary timing results which verify the scalability of the
parallel code.



\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 vectors $\UU$ and $\nn$
$$
H(\UU,\UU;\ \nn) = \EE(\UU) n_1 + \GG(\UU) n_2 .
$$

Then
$$
\HH_{ij}  = \HH(\overline{\UU}_i(m_{ij})),
\overline{\UU}_i(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 on the
\textit{ij}th edge. The function $\overline{\UU}_i(x)$ for $x \in
T_i$ is obtained from the average values, $\UU_i$, of 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{Validation}
\label{sec:validation} The process of validating the \AnuGA{}
application is in its early stages, however initial indications are
encouraging.

As part of the Third International Workshop on Long-wave Runup
Models in 2004 (\url{http://www.cee.cornell.edu/longwave}), four
benchmark problems were specified to allow the comparison of
numerical, analytical and physical models with laboratory and field
data. One of these problems describes a wave tank simulation of the
1993 Okushiri Island tsunami off Hokkaido, Japan \cite{MatH2001}. A
significant feature of this tsunami was a maximum run-up of 32~m
observed at the head of the Monai Valley. This run-up was not
uniform along the coast and is thought to have resulted from a
particular topographic effect. Among other features, simulations of
the Hokkaido tsunami should capture this run-up phenomenon.

\begin{figure}[htbp]
\centerline{\includegraphics[width=4in]{tsunami-fig-3.eps}}
\caption{Comparison of wave tank and \AnuGA{} water stages at gauge
5.}\label{fig:val}
\end{figure}


\begin{figure}[htbp]
\centerline{\includegraphics[width=4in]{tsunami-fig-4.eps}}
\caption{Complex reflection patterns and run-up into Monai Valley
simulated by \AnuGA{} and visualised using our netcdf OSG
viewer.}\label{fig:run}
\end{figure}


The wave tank simulation of the Hokkaido tsunami was used as the
first scenario for validating \AnuGA{}. The dataset provided
bathymetry and topography along with initial water depth and the
wave specifications. The dataset also contained water depth time
series from three wave gauges situated offshore from the simulated
inundation area.

Figure~\ref{fig:val} compares the observed wave tank and modelled
\AnuGA{} water depth (stage height) at one of the gauges. The plots
show good agreement between the two time series, with \AnuGA{}
closely modelling the initial draw down, the wave shoulder and the
subsequent reflections. The discrepancy between modelled and
simulated data in the first 10 seconds is due to the initial
condition in the physical tank not being uniformly zero. Similarly
good comparisons are evident with data from the other two gauges.
Additionally, \AnuGA{} replicates exceptionally well the 32~m Monai
Valley run-up, and demonstrates its occurrence to be due to the
interaction of the tsunami wave with two juxtaposed valleys above
the coastline. The run-up is depicted in Figure~\ref{fig:run}.

This successful replication of the tsunami wave tank simulation on a
complex 3D beach is a positive first step in validating the \AnuGA{}
modelling capability. Subsequent validation will be conducted as
additional datasets become available.



\section{Parallel Implementation}\label{sec:parallel}

To parallelize our algorithm we use a simple domain decomposition
strategy. We suppose we have a global mesh which defines the domain
of our problem. We must first subdivide the global mesh into a set
of non-overlapping submeshes. This partitioning is done using the
\Metis{} partitioning library. We will demonstrate the efficiency of
this library in the following subsections. Once this partitioning
has been done, and the proper communication patterns set, we can
start to evolve our solution. Each timestep consists of
independently and then communicate the appropriate ``ghost'' cell
information once each timestep.

%\begin{figure}[hbtp]
%  \centerline{ \includegraphics[scale = 0.6]{domain.eps}}
%  \caption{The first step of the parallel algorithm is to divide
%  the mesh over the processors.}
%  \label{fig:subpart}
%\end{figure}



\begin{figure}[h!]
  \centerline{ \includegraphics[width=6.5cm]{mermesh.eps}}
  \caption{The global Merimbula mesh}
 \label{fig:mergrid}
\end{figure}

\begin{figure}[h!]
  \centerline{ \includegraphics[width=6.5cm]{mermesh4c.eps}
  \includegraphics[width=6.5cm]{mermesh4a.eps}}

  \centerline{ \includegraphics[width=6.5cm]{mermesh4d.eps}
  \includegraphics[width=6.5cm]{mermesh4b.eps}}
  \caption{The global Merimbula mesh partitioned into 4 submeshes using \Metis.}
 \label{fig:mergrid4}
\end{figure}




\begin{table}[h!]
\caption{4-way and 8-way partition tests of Merimbula mesh showing
the percentage distribution of cells between the
submeshes.}\label{tbl:mer}
\begin{center}
\begin{tabular}{|c|c c c c|}
\hline
CPU & 0 & 1 & 2 & 3\\
\hline
Cells & 2757 & 2713 & 2761 & 2554\\
\% & 25.6\% & 25.2\% & 25.6\% & 23.7\%\ \\
\hline
\end{tabular}

\medskip

\begin{tabular}{|c|c c c c c c c c|}
\hline
CPU & 0 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\
\hline
Cells & 1229 & 1293 & 1352 & 1341 & 1349 & 1401 & 1413 & 1407\\
\% & 11.4\% & 12.0\% & 12.5\% & 12.4\% & 12.5\% & 13.0\% & 13.1\% & 13.1\%\\
\hline
\end{tabular}
\end{center}
\end{table}


\section{Conclusions}
\label{sec:6} \AnuGA{} is a flexible and robust modelling system
that simulates hydrodynamics by solving the shallow water wave
equation in a triangular mesh. It can model the process of wetting
and drying as water enters and leaves an area and is capable of
capturing hydraulic shocks due to the ability of the finite-volume
method to accommodate discontinuities in the solution.

\AnuGA{} can take as input bathymetric and topographic datasets and
simulate the behaviour of riverine flooding, storm surge and
tsunami. Initial validation using wave tank data supports AnuGA's
ability to model complex scenarios. Further validation will be
pursued as additional datasets become available.

\AnuGA{} is already being used to model the behaviour of
hydrodynamic natural hazards. This modelling capability is part of
Geoscience Australia's ongoing research effort to model and
understand the potential impact from natural hazards in order to
reduce their impact on Australian communities.


\paragraph{Acknowledgements:}
The authors are grateful to Belinda Barnes, National Centre for
Epidemiology and Population Health, Australian National University,
and Matt Hayne and Augusto Sanabria, Risk Research Group, Geoscience
Australia, for helpful reviews of a previous version of this paper.
Author Nielsen publish with the permission of the CEO, Geoscience
Australia.

% Preferably provide your bibliography as a separate .bbl file.
% Include URLs, DOIs, Math Review numbers or Zentralblatt numbers in your bibliography
% so we help connect the web of science and ensure maximum visibility
% for your article.
\bibliographystyle{plain}
\bibliography{database1}




\end{document}