Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/numerical_varying_width.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/numerical_varying_width.py	(revision 8886)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/numerical_varying_width.py	(revision 8886)
@@ -0,0 +1,149 @@
+"""Simple water flow example using ANUGA
+
+Water driven up a linear slope and time varying boundary,
+similar to a beach environment
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+from math import cos
+from numpy import zeros, ones, array, interp, polyval, ones_like, zeros_like
+from numpy import where, logical_and
+from time import localtime, strftime, gmtime
+from scipy.interpolate import interp1d
+from anuga.geometry.polygon import inside_polygon, is_inside_triangle
+#from balanced_dev import *
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+#output_dir = 'varying_width'+time
+output_dir = '.'
+output_file = 'varying_width'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1.
+dy = dx
+L = 1500.
+W = 60.
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.,-W/2.))
+
+#domain = anuga.Domain(points, vertices, boundary) 
+domain = Domain(points, vertices, boundary) 
+
+domain.set_name(output_file)                
+domain.set_datadir(output_dir) 
+
+#------------------------------------------------------------------------------
+# Setup Algorithm, either using command line arguments
+# or override manually yourself
+#------------------------------------------------------------------------------
+from anuga.utilities.argparsing import parse_standard_args
+alg, cfl = parse_standard_args()
+domain.set_flow_algorithm(alg)
+domain.set_CFL(cfl)
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+domain.set_quantity('friction', 0.0)
+domain.set_quantity('stage', 12.0)
+XX = array([0.,50.,100.,150.,250.,300.,350.,400.,425.,435.,450.,470.,475.,500.,
+            505.,530.,550.,565.,575.,600.,650.,700.,750.,800.,820.,900.,950.,
+            1000.,1500.])
+ZZ = array([0.,0.,2.5,5.,5.,3.,5.,5.,7.5,8.,9.,9.,9.,9.1,9.,9.,6.,5.5,5.5,5.,
+            4.,3.,3.,2.3,2.,1.2,0.4,0.,0.])
+WW = array([40.,40.,30.,30.,30.,30.,25.,25.,30.,35.,35.,40.,40.,40.,45.,45.,50.,
+            45.,40.,40.,30.,40.,40.,5.,40.,35.,25.,40.,40.])/2.
+
+
+depth = interp1d(XX, ZZ)
+width = interp1d(XX, WW)
+
+def bed_elevation(x,y):
+    z = 25.0*ones_like(x)
+    wid = width(x)
+    dep = depth(x)
+    z = where( logical_and(y < wid, y>-wid), dep, z)
+    return z
+
+domain.set_quantity('elevation', bed_elevation)
+
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+from math import sin, pi, exp
+Br = anuga.Reflective_boundary(domain)      # Solid reflective wall
+#Bt = anuga.Transmissive_boundary(domain)    # Continue all values on boundary 
+#Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+##from anuga.visualiser import RealtimeVisualiser
+##vis = RealtimeVisualiser(domain)
+##vis.render_quantity_height("stage", zScale =h0*500, dynamic=True)
+##vis.colour_height_quantity('stage', (0.0, 0.5, 1.0))
+##vis.start()
+#===============================================================================
+
+
+#---------------------------------------------
+# Find triangle that contains the point Point
+# and print to file
+#---------------------------------------------
+##Point = (0.0, 0.0)
+##for n in range(len(domain.triangles)):
+##    tri = domain.get_vertex_coordinates(n)
+##    if is_inside_triangle(Point,tri):
+##        #print 'Point is within triangle with vertices '+'%s'%tri
+##        n_point = n
+##        break
+##print 'n_point = ',n_point
+##bla
+
+
+
+
+#------------------------------------------------------------------------------
+# Produce a documentation of parameters
+#------------------------------------------------------------------------------
+parameter_file=open('parameters.tex', 'w')
+parameter_file.write('\\begin{verbatim}\n')
+from pprint import pprint
+pprint(domain.get_algorithm_parameters(),parameter_file,indent=4)
+parameter_file.write('\\end{verbatim}\n')
+parameter_file.close()
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 0.1, finaltime = 5.0):
+    #print domain.timestepping_statistics(track_speeds=True)
+    print domain.timestepping_statistics()
+    #vis.update()
+
+
+#test against know data
+    
+#vis.evolveFinished()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/plot_results.py	(revision 8886)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/plot_results.py	(revision 8886)
@@ -0,0 +1,53 @@
+"""
+    Quick plot of the dam break outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from numpy import ones, zeros, ones_like
+
+p_st = util.get_output('varying_width.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[116]
+v2=(p2_st.y==v)
+
+
+#p_dev = util.get_output('dam_break.sww', 0.001)
+#p2_dev=util.get_centroids(p_dev, velocity_extrapolation=True)
+
+#Plot stages
+pyplot.clf()
+pyplot.plot(p2_st.x[v2], p2_st.stage[-1,v2],'b.', label='numerical stage')
+pyplot.plot(p2_st.x[v2], 12.*ones_like(p2_st.x[v2]),'r--', label='analytical stage')
+pyplot.plot(p2_st.x[v2], p2_st.elev[v2],'k-', label='discretised bed')
+pyplot.title('Stage at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.xlim([0,1500])
+pyplot.savefig('stage_plot.png')
+
+
+#Plot xmomentum
+pyplot.clf()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[-1,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], zeros(len(p2_st.x[v2])),'r-', label='analytical')
+pyplot.title('Xmomentum at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.xlim([0,1500])
+pyplot.savefig('xmom_plot.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[-1,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], zeros(len(p2_st.x[v2])),'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.xlim([0,1500])
+pyplot.savefig('xvel_plot.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/produce_report.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/produce_report.py	(revision 8886)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/produce_report.py	(revision 8886)
@@ -0,0 +1,21 @@
+#! /usr/bin/python
+
+# To change this template, choose Tools | Templates
+# and open the template in the editor.
+
+__author__="steve"
+__date__ ="$13/03/2013 4:24:51 PM$"
+
+if __name__ == "__main__":
+    print "Hello World";
+
+import os
+
+
+os.system('pdflatex report.tex')
+os.system('bibtex report')
+os.system('pdflatex report.tex')
+os.system('bibtex report')
+os.system('pdflatex report.tex')
+os.system('bibtex report')
+os.system('pdflatex report.tex')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/produce_results.py	(revision 8886)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/produce_results.py	(revision 8886)
@@ -0,0 +1,19 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from anuga_validation_tests.utilities.fabricate import *
+from anuga_validation_tests.utilities.run_validation import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_steep_island.py')
+    run_validation_script('plot_results.py')
+    run('python', 'produce_report.py')      
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/report.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/report.tex	(revision 8886)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/report.tex	(revision 8886)
@@ -0,0 +1,39 @@
+\documentclass[11pt,a4paper]{article}%{amsart}
+
+\usepackage{graphicx}
+\usepackage{epstopdf}
+\usepackage[section]{placeins} % 'one-shot' command to nicely place figures
+\usepackage{datatool}
+\usepackage{url}
+\DeclareUrlCommand\UScore{\urlstyle{rm}}
+
+
+\input{../../local-defs}
+\input{../../saved_parameters}
+
+%=========================================
+\begin{document} 
+%=========================================
+
+\title{Automated Report on the Performance of \anuga{} on the Steady State of a River at Rest with Varying Width and Topography}
+\maketitle
+%\tableofcontents
+
+%======================
+%\chapter{Dry Dam Break}
+%======================
+
+\input{results}
+
+\section{Parameters used for this simulation}
+
+\input{parameters}
+
+%======================
+% bibliography
+%======================
+
+\bibliographystyle{plain}
+\bibliography{../../bibliography}
+
+\end{document}
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/results.tex	(revision 8886)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/river_at_rest_varying_topo_width/results.tex	(revision 8886)
@@ -0,0 +1,76 @@
+
+\section{River at rest with varying width and topography}
+
+This is a test if the method is well-balanced. Furthermore, we test if the wet/dry interface has been correctly treated for steep river banks. This test is taken from the work of Goutal and Maurel~\cite{GM1997}.
+
+The initial condition is a river at rest with water depth $12.0$. Boundary conditions are solid wall. The width and topography of the river are defined as piecewise linear interpolations of the data presented in Table~\ref{tab:width_topo}.
+\begin{table}
+	\centering
+	\caption{River width and topography (see Goutal and Maurel~\cite{GM1997}).}
+	\begin{tabular}{c|c|c}	\hline
+	$x$&      topography  & width                                                       \\	\hline \hline
+0 &0& 40 \\
+50 &0& 40\\
+100 &2.5& 30\\
+150 &5& 30\\
+250 &5& 30\\
+300 &3& 30\\
+350 &5& 25\\
+400 &5& 25\\
+425 &7.5& 30\\
+435 &8& 35\\
+450 &9& 35\\
+470 &9& 40\\
+475 &9& 40\\
+500 &9.1& 40\\
+505 &9& 45\\
+530 &9& 45\\
+550 &6& 50\\
+565 &5.5& 45\\
+575 &5.5& 40\\
+600 &5& 40\\
+650 &4& 30\\
+700 &3& 40\\
+750 &3& 40\\
+800 &2.3& 5\\
+820 &2& 40\\
+900 &1.2& 35\\
+950 &0.4& 25\\
+1000 &0& 40\\
+1500 &0& 40 \\ \hline
+	\end{tabular}
+	
+	\label{tab:width_topo}
+\end{table} 
+The analytical solution is the river at rest, that is, $w=12.0$ and $u=v=0$.
+
+\subsection{Results}
+
+
+Current version of \anuga{} might not treat the wet/dry interface appropriately. The following three figures show the stage, $x$-momentum, and $x$-velocity respectively, after we run the simulation for some time. We should see excellent agreement between the analytical and numerical solutions if the method is well-balanced and if the wet/dry interface has been correctly treated.
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
