Index: /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/analytical_subcritical.py
===================================================================
--- /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/analytical_subcritical.py	(revision 9228)
+++ /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/analytical_subcritical.py	(revision 9228)
@@ -0,0 +1,44 @@
+"""
+Subcritical flow over a bump.
+Ref1: Houghton & Kasahara, Nonlinear shallow fluid flow over an isolated ridge.
+Comm. Pure and Applied Math. DOI:10.1002/cpa.3160210103
+
+Ref2: Delestre et al, 2012, SWASHES: a compilation of shallow water
+analytic solutions..., Int J Numer Meth Fluids, DOI:10.1002/fld.3741
+
+Sudi Mungkasi, ANU 2012
+"""
+from numpy import zeros, linspace
+from scipy.optimize import fsolve
+from pylab import plot, ylim, show
+from anuga import g
+
+qA = 4.42  # This is the imposed momentum
+hx = 2.0   # This is the water height downstream
+
+def analytic_sol(x):    
+    def elevation(x):
+        z_b = zeros(len(x))
+        for i in range(len(x)):
+            if (8.0 <= x[i] <= 12.0):
+                z_b[i] = 0.2 - 0.05*(x[i]-10.0)**2.0
+            else:
+                z_b[i] = 0.0
+        return z_b
+    z = elevation(x)
+
+    def find_hL(h): #to find the water height at every spatial point
+        return h**3 + (zb-qA**2/(2*g*hx**2)-hx)*h**2 + qA**2/(2*g)
+    h = zeros(len(x))
+    for i in range(len(x)):
+        zb = z[i]
+        h[i] = fsolve(find_hL, 2.0)
+    return h,z
+
+##N = 401
+##L = 25.
+##x = linspace(0.0,L,N)
+##h,z = analytic_sol(x)
+##plot(x,h+z, x,z)
+##ylim([-0.1, 2.1])
+##show()
Index: /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/numerical_subcritical.py
===================================================================
--- /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/numerical_subcritical.py	(revision 9228)
+++ /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/numerical_subcritical.py	(revision 9228)
@@ -0,0 +1,115 @@
+"""
+Simple water flow example using ANUGA
+Subcritical flow over a bump.
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+from anuga import myid, finalize, distribute
+from math import cos
+from numpy import zeros, ones, float
+from time import localtime, strftime, gmtime
+#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 = 'subcritical_'+time
+output_dir = '.'
+output_file = 'subcritical'
+
+args = anuga.get_args()
+alg = args.alg
+verbose = args.verbose
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 0.1
+dy = dx
+L = 25.
+W = 3*dx
+
+if myid == 0:
+    # structured mesh
+    points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, 0.0))
+    
+    #domain = anuga.Domain(points, vertices, boundary) 
+    domain = Domain(points, vertices, boundary) 
+    
+    domain.set_name(output_file)                
+    domain.set_datadir(output_dir) 
+    domain.set_flow_algorithm(alg)
+    
+    #------------------------------------------------------------------------------
+    # Setup initial conditions
+    #------------------------------------------------------------------------------
+    def elevation(x,y):
+        z_b = zeros(len(x))
+        for i in range(len(x)):
+            if (8.0 <= x[i] <= 12.0):
+                z_b[i] = 0.2 - 0.05*(x[i]-10.0)**2.0
+            else:
+                z_b[i] = 0.0
+        return z_b
+    domain.set_quantity('elevation',elevation)
+    domain.set_quantity('friction', 0.0)
+    
+    
+    def height(x,y):
+        return 2.0*ones(len(x))
+    domain.set_quantity('stage', height)
+else:
+    domain = None
+    
+#---------------------------
+# Create Parallel Domain
+#---------------------------    
+domain = distribute(domain)
+
+#-----------------------------------------------------------------------------
+# 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([2., 4.42, 0.]) # Constant boundary values
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bd, 'right': Bd, 'top': Br, 'bottom': Br})
+
+
+
+
+#------------------------------------------------------------------------------
+# Produce a documentation of parameters
+#------------------------------------------------------------------------------
+if myid == 0:
+    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 = 1.0, finaltime = 100.):
+    #print domain.timestepping_statistics(track_speeds=True)
+    if myid == 0 and verbose: print domain.timestepping_statistics()
+
+
+domain.sww_merge(delete_old=True)
+
+
+finalize()
Index: /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/plot_results.py
===================================================================
--- /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/plot_results.py	(revision 9228)
+++ /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/plot_results.py	(revision 9228)
@@ -0,0 +1,55 @@
+"""
+    Quick plot of the subcritical-flow outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_subcritical import *
+from numpy import ones, arange
+
+p_st = util.get_output('subcritical.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+#v2=(p2_st.y==v)
+v#2=(p2_st.y>-1.0)
+v2 = arange(len(p2_st.y))
+
+h,z = analytic_sol(p2_st.x[v2])
+
+#Plot the stages##############################################################
+pyplot.clf()
+pyplot.plot(p2_st.x[v2], p2_st.stage[-1,v2], 'b.-', label='numerical stage') # 0*T/6
+pyplot.plot(p2_st.x[v2], h+z,'r-', label='analytical stage')
+pyplot.plot(p2_st.x[v2], z,'k-', label='bed elevation')
+pyplot.title('Stage at time = %s secs'% p2_st.time[-1])
+##pyplot.ylim(-5.0,5.0)
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png')
+
+
+#Plot the momentums##########################################################
+pyplot.clf()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[-1,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 4.42*ones(len(p2_st.x[v2])),'r-', label='analytical')
+pyplot.title('Xmomentum at time = %s secs'% p2_st.time[-1])
+pyplot.legend(loc='best')
+##pyplot.ylim([1.52,1.54])
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('xmom_plot.png')
+
+
+
+#Plot the velocities#########################################################
+pyplot.clf()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[-1,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 4.42/h,'r-', label='analytical')
+pyplot.title('Xvelocity at time = %s secs'% p2_st.time[-1])
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
Index: /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/produce_results.py
===================================================================
--- /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/produce_results.py	(revision 9228)
+++ /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/produce_results.py	(revision 9228)
@@ -0,0 +1,8 @@
+import anuga
+from anuga.validation_utilities import produce_report
+
+args = anuga.get_args()
+
+produce_report('numerical_subcritical.py', args=args)
+
+
Index: /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/report.tex
===================================================================
--- /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/report.tex	(revision 9228)
+++ /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/report.tex	(revision 9228)
@@ -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{../../reports/local-defs}
+\input{../../reports/saved_parameters}
+
+%=========================================
+\begin{document} 
+%=========================================
+
+\title{Automated Report on the Performance of \anuga{} on the Steady State of Moving Water with an Immersed Bump:\\Subcritical Flow without a Shock}
+\maketitle
+%\tableofcontents
+
+%======================
+%\chapter{Dry Dam Break}
+%======================
+
+\input{results}
+
+\section{Parameters used for this simulation}
+
+\input{parameters}
+
+%======================
+% bibliography
+%======================
+
+\bibliographystyle{plain}
+\bibliography{../../reports/bibliography}
+
+\end{document}
Index: /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/results.tex
===================================================================
--- /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/results.tex	(revision 9228)
+++ /trunk/anuga_core/validation_tests/analytical_exact/supercritical_over_bump/results.tex	(revision 9228)
@@ -0,0 +1,59 @@
+
+\section{Subcritical flow without a shock over a bump}
+
+This is a subcritical flow over a bump. This test is adapted from Goutal and Maurel~\cite{GM1997}. No shock occurs in this scenario. 
+
+Consider a one dimensional domain $[0,25]$ with topography
+\begin{equation}
+z(x)= \left\{ \begin{array}{ll}
+ 0.2-0.05\left(x-10\right)^2& ~\textrm{if}\quad 8 \leq x \leq 12\,,\\
+ 0 & ~\textrm{otherwise}\,,\\
+\end{array} \right.
+\end{equation}
+together with Dirichlet boundary conditions.
+Physically, the boundary conditions mean that there is a source of flow upstream at the point $x=0^{-}$ and at the same time there exists a sink of flow downstream at the point $x=25^{+}$\,.
+
+
+The analytical height is found by solving the Bernoulli equation. The simplified Bernoulli equation is the following cubic equation
+\begin{equation}
+h^3 + \left(z - \frac{q^2}{2 g H^2} - H \right) h^2 + \frac{q^2}{2 g} = 0\,,
+\end{equation}
+where $H$ is the upstream height and $q=uh$ is the discharge or $x$-momentum. When the height $h$ has been found, the velocity is computed as $u=q/h$\,.
+
+\subsection{Results}
+For our test we consider the initial condition
+\begin{equation}
+u(x,y,0)=v(x,y,0)=0\,, \quad
+w(x,y,0)= 0.2\,,
+\end{equation}
+and the Dirichlet boundary conditions at $x=0^{-}$ and $25^{+}$ to be 
+\begin{equation}
+[w,hu,hv]=[2, 4.42, 0]\,.
+\end{equation}
+Representatives of the simulation results are given in the following three figures. Even though we have small discrepancy in the numerical and analytical momenta, these numerical an analytical solutions should agree quite well.
+
+\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
