Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/analytical_MacDonald.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/analytical_MacDonald.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/analytical_MacDonald.py	(revision 8734)
@@ -0,0 +1,89 @@
+"""
+MacDonald's steady flow on a short channel.
+
+Ref1: MacDonald I, Baines MJ, Nichols NK, Samuels PG. Analytic benchmark solutions for open-channel
+flows. Journal of Hydraulic Engineering November 1997; 123(11):1041-1045.
+DOI: 10.1061/(ASCE)0733-9429(1997)123:11(1041).
+
+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, array, sqrt
+from scipy.integrate import quad
+from pylab import plot,show,ylim
+
+
+L     = 100.  # Channel length
+q = q_s = 2.    # Steady discharge
+g     = 9.81  # Accelleration due to gravity
+q_up  = q_s   # Upstream discharge
+#h_down= h_L   # Downstream water height
+
+#Parameters for water height
+a1 = 0.674202
+a2 = 21.7112
+a3 = 14.492
+a4 = 1.4305
+x_s=200./3   # Specified position of a shock
+n  = 0.0328  # Manning's friction coefficient
+
+
+def water_height(x):
+    N = len(x)
+    h = zeros(N)
+    for i in range(N):
+        if x[i] <= x_s:
+            h[i] = (4./g)**(1./3) * (4./3 - x[i]/100.) - 9.*x[i]/1000*(x[i]/100 - 2./3)
+        else:
+            r = x[i]/100. - 2./3.
+            h[i] = (4./g)**(1./3) * (a1*r**4 + a1*r**3 - a2*r**2 + a3*r + a4)
+    return h
+
+##def diff_water_height(x):
+##    N = len(x)
+##    hx = zeros(N)
+##    for i in range(N):
+##        if x[i] <= x_s:
+##            hx[i] = (4./g)**(1./3)*(-0.01) - (9./1000)*(x[i]/100 - 2./3) - (9.*x[i]/1000)*0.01
+##        else:
+##            r  = x[i]/100.-2./3
+##            rx = 0.01
+##            hx[i] = (4./g)**(1./3)*(a1*4.*r**3*rx + a1*3.*r**2*rx - a2*2.*r*rx + a3*rx)
+##    return hx
+
+def bed(x):
+    N = len(x)
+    Z = zeros((N,2))
+    def func1(x):
+        h = (4./g)**(1./3) * (4./3 - x/100.) - 9.*x/1000*(x/100 - 2./3)
+        hx = (4./g)**(1./3)*(-0.01) - (9./1000)*(x/100 - 2./3) - (9.*x/1000)*0.01                
+        return (q**2/(g*h**3) - 1.)*hx - n**2*q*abs(q)/h**(10./3)
+    def func2(x):
+        r = x/100. - 2./3.
+        rx = 0.01
+        h = (4./g)**(1./3) * (a1*r**4 + a1*r**3 - a2*r**2 + a3*r + a4)
+        hx = (4./g)**(1./3)*(a1*4.*r**3*rx + a1*3.*r**2*rx - a2*2.*r*rx + a3*rx)                
+        return (q**2/(g*h**3) - 1.)*hx - n**2*q*abs(q)/h**(10./3)    
+    for i in range(N):
+        if x[i] <= x_s:
+            Z[i] = quad(func1, x[i], x_s)
+            Z[i] += quad(func2, x_s, L)
+        else:           
+            Z[i] = quad(func2, x[i], L)
+    #print 'Integral = ', -Z[:,0], ' with error = ', Z[:,1]
+    return -1.*Z[:,0]
+
+def analytic_sol(x):
+    return water_height(x), bed(x)
+
+##N = 1001
+##X = linspace(0.,L,N)
+##H = water_height(X)
+##Hx= diff_water_height(X)
+##Z = bed(X)
+##
+##plot(X,H+Z, X,Z)
+##show()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/numerical_MacDonald.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/numerical_MacDonald.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/numerical_MacDonald.py	(revision 8734)
@@ -0,0 +1,137 @@
+"""
+Simple water flow example using ANUGA
+Transcritical flow over a bump with a shock
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+from math import cos
+from numpy import zeros, ones, array
+from time import localtime, strftime, gmtime
+from scipy.integrate import quad
+#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 = '.'
+output_file = 'MacDonald'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 0.25
+dy = dx
+L = 100.
+W = 3*dx
+
+# 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) 
+
+#------------------------------------------------------------------------------
+# 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
+#------------------------------------------------------------------------------
+L     = 100.  # Channel length
+q = q_s = 2.    # Steady discharge
+g     = 9.81  # Accelleration due to gravity
+q_up  = q_s   # Upstream discharge
+#h_down= h_L   # Downstream water height
+
+#Parameters for water height
+a1 = 0.674202
+a2 = 21.7112
+a3 = 14.492
+a4 = 1.4305
+x_s=200./3   # Specified position of a shock
+n  = 0.0328  # Manning's friction coefficient
+
+
+def bed(x,y):
+    N = len(x)    
+    Z = zeros((N,2))
+    def func1(x):
+        h = (4./g)**(1./3) * (4./3 - x/100.) - 9.*x/1000*(x/100 - 2./3)
+        hx = (4./g)**(1./3)*(-0.01) - (9./1000)*(x/100 - 2./3) - (9.*x/1000)*0.01                
+        return (q**2/(g*h**3) - 1.)*hx - n**2*q*abs(q)/h**(10./3)
+    def func2(x):
+        r = x/100. - 2./3.
+        rx = 0.01
+        h = (4./g)**(1./3) * (a1*r**4 + a1*r**3 - a2*r**2 + a3*r + a4)
+        hx = (4./g)**(1./3)*(a1*4.*r**3*rx + a1*3.*r**2*rx - a2*2.*r*rx + a3*rx)                
+        return (q**2/(g*h**3) - 1.)*hx - n**2*q*abs(q)/h**(10./3)    
+    for i in range(N):
+        if x[i] <= x_s:
+            Z[i] = quad(func1, x[i], x_s)
+            Z[i] += quad(func2, x_s, L)
+        else:           
+            Z[i] = quad(func2, x[i], L)
+    #print 'Integral = ', -Z[:,0], ' with error = ', Z[:,1]
+    elevation = -1.*Z[:,0]
+    return elevation
+
+domain.set_quantity('stage', 2.87870797)
+domain.set_quantity('elevation', bed)
+domain.set_quantity('friction', n)
+
+
+#-----------------------------------------------------------------------------
+# 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 
+BdL = anuga.Dirichlet_boundary([3.58431872, q_s, 0.]) # Constant boundary values
+BdR = anuga.Dirichlet_boundary([2.87870797, q_s, 0.]) # Constant boundary values
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': BdL, 'right': BdR, '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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 1.0, finaltime = 100.):
+    #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/MacDonald_short_channel/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/plot_results.py	(revision 8734)
@@ -0,0 +1,55 @@
+"""
+    Quick plot for outputs of the transcritical flow with a shock
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_MacDonald import *
+from numpy import ones
+
+p_st = util.get_output('MacDonald.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+h,z = analytic_sol(p2_st.x[v2])
+
+#Plot the stages##############################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[100,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 an instant in time')
+##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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[100,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 2.*ones(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.savefig('xmom_plot.png')
+
+
+
+#Plot the velocities#########################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[100,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 2./h,'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_MacDonald.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/MacDonald_short_channel/results.tex	(revision 8734)
@@ -0,0 +1,35 @@
+
+\section{MacDonald's steady flow: transcritical flow with a shock}
+
+This is a MacDonald's steady flow with a shock.
+
+\subsection{Results}
+
+
+We should see excellent agreement between the analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/analytical_avalanche_dry.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/analytical_avalanche_dry.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/analytical_avalanche_dry.py	(revision 8734)
@@ -0,0 +1,89 @@
+"""
+Analytical solution to debris avalanche involving a dry area
+Ref: Mungkasi and Roberts, ANZIAM J (CTAC2010)
+"""
+
+from numpy import array, zeros, linspace
+from numpy import arctan, sqrt, sin, cos, tan
+
+#Parameters
+g = 9.81           # gravity
+h_0 = 20.0         # depth upstream. Note that the depth downstream is 0.0
+L = 200.0          # length of stream/domain
+
+
+#delt = 0.0
+friction_slope = 0.05             #tan(delta) # NOTE THAT friction_slope must less than bed_slope
+#thet = pi/36.0
+bed_slope = 0.1                   #tan(theta) #0.1      # bottom slope, positive if it is increasing bottom.
+thet = arctan(bed_slope)
+F2 = g*cos(thet)*cos(thet)*friction_slope
+c0 = sqrt(g*h_0)                  # sound speed
+m = -1.0*g*bed_slope + F2         # auxiliary variable
+
+
+# Function for analytical solution
+def analytical_sol(X,t):
+    N = len(X)
+    u = zeros(N)
+    h = zeros(N)
+    w = zeros(N)
+    z = zeros(N)
+    mom = zeros(N)
+    if abs(t)>1e-9:
+        for i in range(N):
+            # Calculate Analytical Solution at time t > 0
+            if X[i] <= -2.0*c0*t + 0.5*m*t**2.0:
+                u[i] = 0.0
+                h[i] = 0.0
+            elif X[i] <= c0*t + 0.5*m*t**2.0:
+                u[i] = 2.0/3.0 * (X[i]/t - c0 + m*t)
+                h[i] = 1.0/(9.0*g) * (X[i]/t + 2.0*c0 - 0.5*m*t)**2.0
+            else:
+                u[i] = m*t
+                h[i] = h_0
+            z[i] = bed_slope*X[i]
+            w[i] = h[i] + z[i]
+            mom[i] = u[i]*h[i]
+    else:
+        for i in range(N):
+            # Calculate Analytical Solution at time t > 0
+            if X[i] <= 0.0:
+                u[i] = 0.0
+                h[i] = 0.0
+            else:
+                u[i] = 0.0
+                h[i] = h_0
+            z[i] = bed_slope*X[i]
+            w[i] = h[i] + z[i]
+            mom[i] = u[i]*h[i]    
+    return u,h,w,z,mom
+
+"""
+# Run a simulation
+n = 800            # number of cells
+cell_len = L/n     # length of each cell
+points = zeros(n+1)
+for i in range (n+1):
+    points[i] = i*cell_len - 0.5*L
+X = points
+t = 2.0
+#The following is for plotting the result.
+Uv,Hv,Wv,Zv,Mv = analytical_sol(X,t)
+from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+hold(False)
+clf()
+
+plot1 = subplot(211)
+plot(X,Wv,'b-', X,Zv,'k-')
+ylabel('Stage')
+#plot1.set_ylim([-50.0,75.0])#([-1.0,21.0])
+legend(('analytical solution', 'bed elevation'), 'upper left', shadow=False)
+
+plot2 = subplot(212)
+plot(X,Mv,'b-')
+xlabel('Position')
+ylabel('Momentum')
+#plot2.set_ylim([-200.0,10.0])#([-90.0,10.0])   
+show()
+"""
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/numerical_avalanche_dry.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/numerical_avalanche_dry.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/numerical_avalanche_dry.py	(revision 8734)
@@ -0,0 +1,172 @@
+"""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
+from time import localtime, strftime, gmtime
+from numpy import sin, cos, tan, arctan
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+#output_dir = 'avalanche_'+time
+output_dir = '.'
+output_file = 'avalanche'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+g = 9.81           # gravity
+h_0 = 20.0         # depth upstream. Note that the depth downstream is 0.0
+L = 200.0          # length of stream/domain
+dx = 0.5
+dy = dx
+W = 3*dx
+
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (-L/2.0, -W/2.0))
+
+#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
+#------------------------------------------------------------------------------
+
+# No Mannings friction, but we introduce Coulomb friction below.
+domain.set_quantity('friction', 0.0)
+class Coulomb_friction:
+    
+    def __init__(self,
+                 friction_slope=0.05,
+                 bed_slope=0.1):
+        
+        self.friction_slope = friction_slope   #tan(delta) # NOTE THAT friction_slope must less than bed_slope
+        self.bed_slope = bed_slope             #tan(theta) # bottom slope, positive if it is increasing bottom.
+
+        thet = arctan(bed_slope)
+        self.F = g*cos(thet)*cos(thet)*friction_slope
+        self.m = -1.0*g*bed_slope + self.F
+        
+    def __call__(self, domain):        
+
+        w = domain.quantities['stage'].centroid_values
+        z = domain.quantities['elevation'].centroid_values
+        h = w-z
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        xmom_update[:] = xmom_update + self.F*h
+
+
+bed_slope = 0.1
+friction_slope = 0.05
+Coulomb_forcing_term = Coulomb_friction(friction_slope, bed_slope)
+domain.forcing_terms.append(Coulomb_forcing_term)
+
+def stage(X,Y):
+    N = len(X)
+    w = zeros(N)
+    for i in range(N):
+        if X[i]<=0.0:
+            w[i] = bed_slope*X[i]
+        else:
+            w[i] = bed_slope*X[i] + h_0
+    return w
+domain.set_quantity('stage', stage)
+
+def elevation(X,Y):
+    N = len(X)
+    y=zeros(N)
+    for i in range(N):
+        y[i] = bed_slope*X[i]
+    return y
+domain.set_quantity('elevation',elevation)
+
+def f_right(t):
+    z_r = bed_slope*(0.5*L)
+    h_r = h_0 #+ bed_slope*cell_len
+    w_r = z_r + h_r
+    u_r = Coulomb_forcing_term.m*t
+    #['stage', 'xmomentum', 'ymomentum']
+    return [w_r,  u_r*h_r,  0.0]
+
+
+
+
+#-----------------------------------------------------------------------------
+# 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
+BTime = anuga.Time_boundary(domain,f_right)
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bt, 'right': BTime, '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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# 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 = 3.):
+    #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/avalanche_dry/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/plot_results.py	(revision 8734)
@@ -0,0 +1,66 @@
+"""
+    Quick plot of the dam break outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_avalanche_dry import *
+
+p_st = util.get_output('avalanche.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+u0,h0,w0,z0,p0 = analytical_sol(p2_st.x[v2], p2_st.time[0])
+u10,h10,w10,z10,p10 = analytical_sol(p2_st.x[v2], p2_st.time[10])
+u30,h30,w30,z30,p30 = analytical_sol(p2_st.x[v2], p2_st.time[30])
+
+#Plot stages
+pyplot.clf()
+#pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[0,v2],'b.-', label='numerical stage')
+pyplot.plot(p2_st.x[v2], p2_st.stage[10,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], p2_st.stage[30,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], w0,'r-', label='analytical stage')
+pyplot.plot(p2_st.x[v2], w10,'r-')
+pyplot.plot(p2_st.x[v2], w30,'r-')
+pyplot.plot(p2_st.x[v2], p2_st.elev[v2],'k-', label='bed elevation')
+pyplot.title('Stage at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png', dpi=600)
+#pyplot.show()
+
+
+#Plot xmomentum
+pyplot.clf()
+#pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[0,v2], 'b.-', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[10,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[30,v2],'b.-')
+pyplot.plot(p2_st.x[v2], u0*h0,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u10*h10,'r-')
+pyplot.plot(p2_st.x[v2], u30*h30,'r-')
+pyplot.title('Xmomentum at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('xmom_plot.png')
+
+
+#Plot velocities
+pyplot.clf()
+#pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[0,v2], 'b.-', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[10,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[30,v2],'b.-')
+pyplot.plot(p2_st.x[v2], u0,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u10,'r-')
+pyplot.plot(p2_st.x[v2], u30,'r-')
+pyplot.title('Xvelocity at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_avalanche_dry.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/report.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/report.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/report.tex	(revision 8734)
@@ -0,0 +1,39 @@
+\documentclass[11pt,a4paper]{article}
+
+\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 Dry Debris Avalanche Problem}
+\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/avalanche_dry/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_dry/results.tex	(revision 8734)
@@ -0,0 +1,67 @@
+
+\section{Avalanche involving a dry area}
+
+An avalanche problem involving a dry area is solved using shallow water approach. This problem is very similar to the dry dam break, but it is on a sloping topography. The debris could be snow, sand, or even rock. The simulation should show a rarefaction and wetting process, just like the dry dam break problem. The analytical solution of this problem was derived by Mungkasi and Roberts~\cite{MR2011DA}. This shallow water approach to solve debris avalanche problems was also implemented by a number of researchers, such as Mangeney et al.~\cite{MHR2000} and Naaim et al.~\cite{NVC1997}. 
+
+The initial condition is
+\begin{equation} \label{eq:dap_init}
+u(x,0)=0, ~~v(x,y)=0, ~~\textrm{and}~~
+h(x,0) = \left\{ \begin{array}{ll}
+h_1 & \textrm{if $x < 0$}\\
+0 & \textrm{if $x > 0$}\\
+\end{array} \right.
+\end{equation}
+where $h_1>0$. The topography is a flat bed with positive slope.
+
+The analytical solution~\cite{MR2011DA} at time $t>0$ is
+\begin{equation} 
+h(x) = \left\{ \begin{array}{ll}
+0 & \textrm{if $x \leq -2 c_0 t + \frac12 mt^2$}\\
+h_R=\frac{1}{9g} \left( \frac{x}{t} + 2c_0 - \frac12 mt \right)^2 & \textrm{if $-2 c_0 t + \frac12 mt^2 \leq x \leq c_0 t + \frac12 mt^2$}\\
+h_0 & \textrm{if $x \geq c_0 t + \frac12 mt^2$}\\
+\end{array} \right.
+\end{equation}
+which is the free surface and
+\begin{equation} 
+u(x) = \left\{ \begin{array}{ll}
+0 & \textrm{if $x \leq -2 c_0 t + \frac12 mt^2$}\\
+u_R=\frac23 \left( \frac{x}{t} - c_0 + mt \right) & \textrm{if $-2 c_0 t + \frac12 mt^2 \leq x \leq c_0 t + \frac12 mt^2$}\\
+mt & \textrm{if $x \geq c_0 t + \frac12 mt^2$}\\
+\end{array} \right.
+\end{equation}
+which is the velocity. Here $m=-g\tan{\theta}+F$, where $\tan{\theta}$ is the slope of the topography. Variable $F$ is the Coulomb-type friction given by 
+\begin{equation}
+F=g \cos^2{\theta} \tan{\delta},
+\end{equation}
+in which $\tan{\delta}$ is a given value of friction slope such that $\tan{\delta} \leq \tan{\theta}$.
+
+
+\subsection{Results}
+
+For our test, we consider $h_0=20$ in (\ref{eq:dap_init}).
+The following figures show the stage, $x$-momentum, and $x$-velocity at several instants of time. We should see excellent agreement between the analytical and numerical solutions. The wet/dry interface is difficult to resolve and it usually produces large errors, similar to the dry dam break problem.
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/analytical_avalanche_wet.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/analytical_avalanche_wet.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/analytical_avalanche_wet.py	(revision 8734)
@@ -0,0 +1,115 @@
+"""
+Analytical solution to debris avalanche involving a shock wave
+Ref: Mungkasi and Roberts, Pure and Applied Geophysics 2012
+"""
+
+from numpy import array, zeros, linspace
+from numpy import arctan, sqrt, sin, cos, tan
+
+#Parameters
+g = 9.81           # gravity
+h_0 = 20.0         # depth upstream
+h_1 = 10.0         # depth downstream
+L = 200.0          # length of stream/domain
+
+
+#delt = 0.0
+friction_slope = 0.05             #tan(delta) # NOTE THAT friction_slope must less than bed_slope
+#thet = pi/36.0
+bed_slope = 0.1                   #tan(theta) #0.1      # bottom slope, positive if it is increasing bottom.
+thet = arctan(bed_slope)
+F2 = g*cos(thet)*cos(thet)*friction_slope
+c0 = sqrt(g*h_0)                  # sound speed
+m = -1.0*g*bed_slope + F2         # auxiliary variable
+
+
+# Function for analytical solution
+def analytical_sol(X,t):
+    #The following is for computing the shock velocity
+    c0 = sqrt(g*h_0) #right celerity
+    c1 = sqrt(g*h_1) #left celerity
+    Smin=-21.0
+    Smax=20.0
+    for i in range(100):
+        S=(Smin+Smax)/2.0       
+        u2=S-c1*c1/4.0/S*(1.0+sqrt(1.0+8.0*S*S/c1/c1))
+        c2=c1*sqrt(0.5*(sqrt(1.0+8.0*S*S/c1/c1)-1.0))
+        func = -2.0*c0 -u2 +2.0*c2
+        if (func > 0.0):
+            Smin=S
+        else:
+            Smax=S
+    if( abs(S) > 99.0):
+        print 'no convergence'
+    #h2= h_1/(1.0-u2/S)
+    h2 = c2**2/9.81
+    #print "S=",S
+    #print "c2=",c2
+    #print "u2=",u2
+    
+    N = len(X)
+    u = zeros(N)
+    h = zeros(N)
+    w = zeros(N)
+    z = zeros(N)
+    mom = zeros(N)
+    if abs(t)>1e-9:
+        for i in range(N):
+            # Calculate Analytical Solution at time t > 0
+            if X[i]-0.5*m*t**2 <= S*t :
+                u[i] = m*t 
+                h[i] = h_1
+            elif X[i]-0.5*m*t**2 <= (u2+c2)*t :
+                u[i] = u2+m*t
+                h[i] = h2
+            elif X[i]-0.5*m*t**2 < c0*t :
+                u[i] = 2.0/3.0*((X[i]-0.5*m*t**2)/t - sqrt(g*h_0))  +m*t
+                h[i] = 4.0/(9.0*g)*((X[i]-0.5*m*t**2)/(2.0*t) + sqrt(g*h_0))**2
+            else:
+                u[i] = m*t
+                h[i] = h_0 
+            z[i] = bed_slope*X[i]
+            w[i] = h[i] + z[i]
+            mom[i] = u[i]*h[i]
+    else:
+        for i in range(N):
+            # Calculate Analytical Solution at time t > 0
+            if X[i] <= 0.0:
+                u[i] = 0.0
+                h[i] = h_1
+            else:
+                u[i] = 0.0
+                h[i] = h_0
+            z[i] = bed_slope*X[i]
+            w[i] = h[i] + z[i]
+            mom[i] = u[i]*h[i]    
+    return u,h,w,z,mom
+
+"""
+# Run a simulation
+n = 800            # number of cells
+cell_len = L/n     # length of each cell
+points = zeros(n+1)
+for i in range (n+1):
+    points[i] = i*cell_len - 0.5*L
+X = points
+t = 4.0
+#The following is for plotting the result.
+Uv,Hv,Wv,Zv,Mv = analytical_sol(X,t)
+from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+hold(False)
+clf()
+
+plot1 = subplot(211)
+plot(X,Wv,'b-', X,Zv,'k-')
+ylabel('Stage')
+#plot1.set_ylim([-50.0,75.0])#([-1.0,21.0])
+legend(('analytical solution', 'bed elevation'), 'upper left', shadow=False)
+
+plot2 = subplot(212)
+plot(X,Mv,'b-')
+xlabel('Position')
+ylabel('Momentum')
+#plot2.set_ylim([-200.0,10.0])#([-90.0,10.0])   
+show()
+"""
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/numerical_avalanche_wet.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/numerical_avalanche_wet.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/numerical_avalanche_wet.py	(revision 8734)
@@ -0,0 +1,182 @@
+"""
+Numerical solution to debris avalanche involving a shock wave
+Ref: Mungkasi and Roberts, Pure and Applied Geophysics 2012
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+from math import cos
+from numpy import zeros
+from time import localtime, strftime, gmtime
+from numpy import sin, cos, tan, arctan
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+#output_dir = 'avalanche_'+time
+output_dir = '.'
+output_file = 'avalanche'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+g = 9.81           # gravity
+h_0 = 20.0         # depth upstream
+h_1 = 10.0         # depth downstream
+L = 200.0          # length of stream/domain
+dx = 0.5
+dy = dx
+W = 3*dx
+
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (-L/2.0, -W/2.0))
+
+#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
+#------------------------------------------------------------------------------
+
+# No Mannings friction, but we introduce Coulomb friction below.
+domain.set_quantity('friction', 0.0)
+class Coulomb_friction:
+    
+    def __init__(self,
+                 friction_slope=0.05,
+                 bed_slope=0.1):
+        
+        self.friction_slope = friction_slope   #tan(delta) # NOTE THAT friction_slope must less than bed_slope
+        self.bed_slope = bed_slope             #tan(theta) # bottom slope, positive if it is increasing bottom.
+
+        thet = arctan(bed_slope)
+        self.F = g*cos(thet)*cos(thet)*friction_slope
+        self.m = -1.0*g*bed_slope + self.F
+        
+    def __call__(self, domain):        
+
+        w = domain.quantities['stage'].centroid_values
+        z = domain.quantities['elevation'].centroid_values
+        h = w-z
+
+        xmom_update = domain.quantities['xmomentum'].explicit_update
+        ymom_update = domain.quantities['ymomentum'].explicit_update
+
+        xmom_update[:] = xmom_update + self.F*h
+
+
+bed_slope = 0.1
+friction_slope = 0.05
+Coulomb_forcing_term = Coulomb_friction(friction_slope, bed_slope)
+domain.forcing_terms.append(Coulomb_forcing_term)
+
+def stage(X,Y):
+    N = len(X)
+    w = zeros(N)
+    for i in range(N):
+        if X[i]<=0.0:
+            w[i] = bed_slope*X[i] + h_1
+        else:
+            w[i] = bed_slope*X[i] + h_0
+    return w
+domain.set_quantity('stage', stage)
+
+def elevation(X,Y):
+    N = len(X)
+    y=zeros(N)
+    for i in range(N):
+        y[i] = bed_slope*X[i]
+    return y
+domain.set_quantity('elevation',elevation)
+
+def f_right(t):
+    z_r = bed_slope*(0.5*L)
+    h_r = h_0 #+ bed_slope*cell_len
+    w_r = z_r + h_r
+    u_r = Coulomb_forcing_term.m*t
+    #['stage', 'xmomentum', 'ymomentum']
+    return [w_r,  u_r*h_r,  0.0]
+
+def f_left(t):
+    z_l = bed_slope*(-0.5*L)
+    h_l = h_1 #+ bed_slope*cell_len
+    w_l = z_l + h_l
+    u_l = Coulomb_forcing_term.m*t
+    #['stage', 'xmomentum', 'ymomentum']
+    return [w_l,  u_l*h_l,  0.0]
+
+
+
+
+#-----------------------------------------------------------------------------
+# 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
+BTimeR = anuga.Time_boundary(domain,f_right)
+BTimeL = anuga.Time_boundary(domain,f_left)
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': BTimeL, 'right': BTimeR, '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()
+#===============================================================================
+
+
+
+#------------------------------------------------------------------------------
+# 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 = 4.):
+    #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/avalanche_wet/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/plot_results.py	(revision 8734)
@@ -0,0 +1,66 @@
+"""
+    Quick plot of the avalanche outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_avalanche_wet import *
+
+p_st = util.get_output('avalanche.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[20]
+v2=(p2_st.y==v)
+
+u0,h0,w0,z0,p0 = analytical_sol(p2_st.x[v2], p2_st.time[0])
+u20,h20,w20,z20,p20 = analytical_sol(p2_st.x[v2], p2_st.time[20])
+u40,h40,w40,z40,p40 = analytical_sol(p2_st.x[v2], p2_st.time[40])
+
+#Plot stages
+pyplot.clf()
+#pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[0,v2],'b.-', label='numerical stage')
+pyplot.plot(p2_st.x[v2], p2_st.stage[20,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], p2_st.stage[40,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], w0,'r-', label='analytical stage')
+pyplot.plot(p2_st.x[v2], w20,'r-')
+pyplot.plot(p2_st.x[v2], w40,'r-')
+pyplot.plot(p2_st.x[v2], p2_st.elev[v2],'k-', label='bed elevation')
+pyplot.title('Stage at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png', dpi=600)
+#pyplot.show()
+
+
+#Plot xmomentum
+pyplot.clf()
+#pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[0,v2], 'b.-', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[20,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[40,v2],'b.-')
+pyplot.plot(p2_st.x[v2], u0*h0,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u20*h20,'r-')
+pyplot.plot(p2_st.x[v2], u40*h40,'r-')
+pyplot.title('Xmomentum at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('xmom_plot.png')
+
+
+#Plot velocities
+pyplot.clf()
+#pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[0,v2], 'b.-', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[20,v2], 'b.-')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[40,v2],'b.-')
+pyplot.plot(p2_st.x[v2], u0,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u20,'r-')
+pyplot.plot(p2_st.x[v2], u40,'r-')
+pyplot.title('Xvelocity at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_avalanche_wet.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/report.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/report.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/report.tex	(revision 8734)
@@ -0,0 +1,38 @@
+\documentclass[11pt,a4paper]{article}
+
+\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 Avalanche involving a Shock}
+\maketitle
+%\tableofcontents
+
+%======================
+%\chapter{Wet 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/avalanche_wet/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/avalanche_wet/results.tex	(revision 8734)
@@ -0,0 +1,85 @@
+
+\section{Avalanche involving a shock wave}
+
+We consider an avalanche involving a shock wave.
+This problem is similar to dam break on wet areas, and so, it involves a shock. We consider a flat topography with positive slope. Shallow water approach is used to solve the problem. The analytical solution of this problem was derived by Mungkasi and Roberts~\cite{MR2012PAAG}. This shallow water approach was also implemented by a number of researchers, such as Mangeney et al.~\cite{MHR2000} and Naaim et al.~\cite{NVC1997}. The simulation should show a rarefaction and a shock. 
+
+
+
+The initial condition is
+\begin{equation} \label{eq:dap_init}
+u(x,0)=0, ~~v(x,y)=0, ~~\textrm{and}~~
+h(x,0) = \left\{ \begin{array}{ll}
+h_1 & \textrm{if $x < 0$}\\
+h_0 & \textrm{if $x > 0$}\\
+\end{array} \right.
+\end{equation}
+where $h_0>h_1>0$. The topography is a flat bed with positive slope. Note that when $h_1=0$, the problem becomes avalanche involving a dry area~\cite{MR2011DA}.
+
+The analytical solution~\cite{MR2012PAAG} at time $t>0$ is
+\begin{equation}
+h(x,t) = \left\{ \begin{array}{ll}
+h_1 & \textrm{if $x < \sigma t + \frac12 mt^2$}\,,\\
+h_2 & \textrm{if $ \sigma t + \frac12 mt^2 \leq x < (u_2+c_2)t + \frac12 mt^2$}\,,\\
+\frac{1}{9g} \left( \frac{x}{t} + 2c_0 - \frac12 mt \right)^2 & \textrm{if $(u_2+c_2)t + \frac12 mt^2 \leq x < c_0 t + \frac12 mt^2$}\,,\\
+h_0 & \textrm{if $x \geq c_0 t + \frac12 mt^2$}\,, \\
+\end{array} \right.
+\end{equation}
+and
+\begin{equation}
+u(x,t) = \left\{ \begin{array}{ll}
+mt & \textrm{if $x < \sigma t + \frac12 mt^2$}\,,\\
+u_2+mt & \textrm{if $ \sigma t + \frac12 mt^2 \leq x < (u_2+c_2)t + \frac12 mt^2$}\,,\\
+\frac23 \left( \frac{x}{t} - c_0 + mt \right) & \textrm{if $(u_2+c_2)t + \frac12 mt^2 \leq x < c_0 t + \frac12 mt^2$}\,,\\
+mt & \textrm{if $x \geq c_0 t + \frac12 mt^2$}\,, \\ \end{array} \right.
+\end{equation}
+for time $t>0$\,. Here $u_2$\,, $c_2$\,, and $\sigma$ are the solutions of the three simultaneous equations 
+\begin{equation}
+u_2  = \sigma - \frac{c_1^2}{4\sigma} \left(1+\sqrt{1+8\left( \frac{\sigma}{c_1} \right)^2}     \right),
+\end{equation}
+\begin{equation}
+c_2 = c_1 \sqrt{\frac12 \left(\sqrt{1+  8 \left(\frac{\sigma}{c_1}\right)^2}-1  \right)},
+\end{equation}
+and 
+\begin{equation}
+-2c_0 = u_2 -2c_2.
+\end{equation} 
+The value of $h_2$ is calculated using relation $c_2=\sqrt{gh_2}$\,.
+Here $m=-g\tan{\theta}+F$, where $\tan{\theta}$ is the slope of the topography. Variable $F$ is the Coulomb-type friction given by 
+\begin{equation}
+F=g \cos^2{\theta} \tan{\delta},
+\end{equation}
+in which $\tan{\delta}$ is a given value of friction slope such that $\tan{\delta} \leq \tan{\theta}$. 
+
+
+
+\subsection{Results}
+
+For our test, we consider $h_0=20$ and $h_1=10$ in (\ref{eq:dap_init}).
+The following figures show the stage, $x$-momentum, and $x$-velocity at several instants of time. We should see excellent agreement between the analytical and numerical solutions. 
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/analytical_carrier_greenspan.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/analytical_carrier_greenspan.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/analytical_carrier_greenspan.py	(revision 8734)
@@ -0,0 +1,127 @@
+from scipy import sin, cos, sqrt, linspace, pi, dot
+from numpy import zeros, array
+from scipy.special import jn
+from scipy.optimize import fsolve
+from gaussPivot import *
+
+
+#Parameters
+g = 9.81
+
+def j0(x):
+    return jn(0.0, x)
+
+def j1(x):
+    return jn(1.0, x)
+
+def analytic_cg(points, t=0.0, h0=5e2, L=5e4, a=1.0, Tp=900.0):
+    #
+    #Exact solution of the Carrier-Greenspan periodic solution.
+    #Ref1: Carrier and Greenspan, J. Fluid Mech., 1958
+    #Ref2: Mungkasi and Roberts, Int. J. Numer. Meth. Fluids, 2012
+    #Here, points are the spatial evaluating points, t is time variable, h0 is the water depth at Origin
+    #L is the horizontal length considered, a is the amplitude of perturbation at Origin
+    #Tp is the period of oscillation.
+    #
+    #Analytical computations#################################################################
+    def shore(t):
+        def g(u):
+            return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u))    
+        u = fsolve(g,0.0)
+        xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u))
+        position = 1.0 + xi
+        return position, u
+
+    def func(q): #Here q=(w, u)
+        f = zeros(2)
+        f[0] = q[0] + 0.5*q[1]**2.0 - A*j0(4.0*pi/T*(1.0+q[0]-x)**0.5)*cos(2.0*pi/T*(tim+q[1]))
+        f[1] = q[1] + A*j1(4.0*pi/T*(1.0+q[0]-x)**0.5)*sin(2.0*pi/T*(tim+q[1]))/(1+q[0]-x)**0.5
+        return f  
+
+    def newtonRaphson2(f,q,tol=1.0e-15): ##1.0e-9 may be too large
+        for i in range(50):                                                                                                                                    
+            h = 1.0e-4      ##1.0e-4 may be too large.
+            n = len(q)
+            jac = zeros((n,n))
+            if 1.0+q[0]-x<0.0:
+                print "Problem occurs i=",i
+                print "PROBLEM OCCURS.......... 1.0+q[0]-x=",1.0+q[0]-x
+                q[0] = x-1.0 + 0.0001
+                print "NOW problem is fixed as  1.0+q[0]-x=",1.0+q[0]-x
+            f0 = f(q)
+            for i in range(n):
+                temp = q[i]
+                q[i] = temp + h
+                f1 = f(q)
+                q[i] = temp
+                jac[:,i] = (f1 - f0)/h
+            if sqrt(dot(f0,f0)/len(q)) < tol: return q
+            dq = gaussPivot(jac,-f0)
+            q = q + dq
+            if sqrt(dot(dq,dq)) < tol*max(max(abs(q)),1.0): return q
+        print 'Too many iterations'    
+    ##################################################################################
+    N = len(points)
+    
+    eps = a/h0           # Dimensionless amplitude
+    T = Tp*sqrt(g*h0)/L  # Dimensionless period
+    tim = t*sqrt(g*h0)/L # Dimensionless time
+    A = eps/j0(4.0*pi/T) # Dimensionless max horizontal displacement of shore
+    
+    W = zeros(N) # Dimensional stage
+    P = zeros(N) # Dimensional momentum
+    Z = zeros(N) # Dimensional bed
+    H = zeros(N) # Dimensional height
+    U = zeros(N) # Dimensional velocity
+    
+    q = zeros(2)
+    pos_shore, vel_shore = shore(tim) # Dimensionless shore
+    #q[0] = pos_shore
+    #q[1] = vel_shore
+
+    for m in range(N):
+        i = N - m-1
+        X = points[i]            # Dimensional
+        if X >= pos_shore*L:     # Dimensional
+            W[i] = (h0/L)*X - h0 # Dimensional
+            U[i] = 0.0           # Dimensional           
+            Z[i] = (h0/L)*X - h0 # Dimensional
+            H[i] = 0.0           # Dimensional
+            P[i] = 0.0           # Dimensional
+        else:
+            x = X/L              # Dimensionless
+            #q = fsolve(func,q)  # Dimensionless (too sensitive and not good)
+            q = newtonRaphson2(func,q)  # Dimensionless (better than fsolve)           
+            W[i] = q[0]*h0      # Dimensional
+            U[i] = q[1]         # It needs dimensionalisation
+            Z[i] = (h0/L)*X - h0# Dimensional
+            H[i] = W[i] - Z[i]  # Dimensional
+            P[i] = H[i]*U[i]    # It needs dimensionalisation
+    U = U*sqrt(g*h0) #This is the dimensionalisation
+    P = P*sqrt(g*h0) #This is the dimensionalisation
+    return W, P, Z, H, U
+
+
+if __name__ == "__main__":
+    points = linspace(0.0, 55000.0, 1000)
+    W, P, Z, H, U = analytic_cg(points,t=300.0, h0=5e2, L=5e4, a=1.0, Tp=900.0)
+
+    from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot,ion
+    hold(False)
+    clf()
+    plot1 = subplot(311)
+    plot(points/1e+4,W,  points/1e+4,Z)
+    ylabel('Stage')
+
+    plot2 = subplot(312)
+    plot(points/1e+4,P,'b-')
+    ylabel('Momentum')      
+
+    plot3 = subplot(313)
+    plot(points/1e+4,U,'b-')
+    xlabel('Position / 10,000')
+    ylabel('Velocity')
+    legend(('Analytical Solution','Numerical Solution'),
+           'upper center', shadow=False)
+
+    show()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/bisect_function.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/bisect_function.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/bisect_function.py	(revision 8734)
@@ -0,0 +1,32 @@
+## module bisect
+''' root = bisect(f,x1,x2,switch=0,tol=1.0e-9).
+    Finds a root of f(x) = 0 by bisection.
+    The root must be bracketed in (x1,x2).
+    Setting switch = 1 returns root = None if
+    f(x) increases as a result of a bisection.
+    
+    Taken from the book Numerical Methods in Engineering with Python
+    by J. Kiusalaas
+'''    
+from math import log,ceil
+import error
+
+def bisect(f,x1,x2,switch=0,epsilon=1.0e-15):#9):
+    f1 = f(x1)
+    if f1 == 0.0: return x1
+    f2 = f(x2)
+    if f2 == 0.0: return x2
+    if f1*f2 > 0.0: error.err('Root is not bracketed')
+    n = ceil(log(abs(x2 - x1)/epsilon)/log(2.0))
+    n = int(n)
+    for i in range(n):
+        x3 = 0.5*(x1 + x2); f3 = f(x3)
+        if (switch == 1) and (abs(f3) >abs(f1)) \
+                         and (abs(f3) > abs(f2)):
+            return None   
+        if f3 == 0.0: return x3
+        if f2*f3 < 0.0:
+            x1 = x3; f1 = f3
+        else:
+            x2 = x3; f2 = f3
+    return (x1 + x2)/2.0
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/error.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/error.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/error.py	(revision 8734)
@@ -0,0 +1,12 @@
+## module error
+''' err(string).
+    Prints 'string' and terminates program.
+
+    Taken from the book Numerical Methods in Engineering with Python
+    by J. Kiusalaas    
+'''    
+import sys
+def err(string):
+    print string
+    raw_input('Press return to exit')
+    sys.exit()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/gaussPivot.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/gaussPivot.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/gaussPivot.py	(revision 8734)
@@ -0,0 +1,46 @@
+## module gaussPivot
+''' x = gaussPivot(a,b,tol=1.0e-9).
+    Solves [a]{x} = {b} by Gauss elimination with
+    scaled row pivoting
+
+    Taken from the book Numerical Methods in Engineering with Python
+    by J. Kiusalaas    
+'''    
+from numpy import *
+import swap 
+import error
+
+def gaussPivot(a,b,tol=1.0e-12):
+    n = len(b)
+    
+  # Set up scale factors
+    s = zeros(n)
+    for i in range(n):
+        s[i] = max(abs(a[i,:]))
+            
+    for k in range(0,n-1):
+        
+      # Row interchange, if needed
+        p = int(argmax(abs(a[k:n,k])/s[k:n])) + k
+        if abs(a[p,k]) < tol: error.err('Matrix is singular')
+        if p != k:
+            swap.swapRows(b,k,p)
+            swap.swapRows(s,k,p)
+            swap.swapRows(a,k,p)
+            
+      # Elimination
+        for i in range(k+1,n):
+            if a[i,k] != 0.0:
+                lam = a[i,k]/a[k,k]
+                a[i,k+1:n] = a [i,k+1:n] - lam*a[k,k+1:n]
+                b[i] = b[i] - lam*b[k]
+    if abs(a[n-1,n-1]) < tol: error.err('Matrix is singular')
+                   
+  # Back substitution
+    for k in range(n-1,-1,-1):
+        b[k] = (b[k] - dot(a[k,k+1:n],b[k+1:n]))/a[k,k]
+    return b
+
+
+
+        
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/numerical_carrier_greenspan.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/numerical_carrier_greenspan.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/numerical_carrier_greenspan.py	(revision 8734)
@@ -0,0 +1,172 @@
+"""
+Periodic water flows  using ANUGA,
+where water driven up a linear sloping beach and time varying boundary.
+Ref1: Carrier and Greenspan, Journal of Fluid Mechanics, 1958
+Ref2: Mungkasi and Roberts, Int. J. Numerical Methods in Fluids, 2012
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+from math import cos
+from numpy import zeros, array, float
+from time import localtime, strftime, gmtime
+from scipy.optimize import fsolve
+from math import sin, pi, exp, sqrt
+from scipy.special import jn
+#from run_carrier_greenspan_analytic import *
+#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 = 'carrier_greenspan_'+time
+output_dir = '.'
+output_file = 'carrier_greenspan'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+#DIMENSIONAL PARAMETERS
+dx  = 100.
+dy  = dx
+L   = 5e4         # Length of channel (m)
+W   = 5*dx        # Width of channel (m)       
+h0 = 5e2          # Height at origin when the water is still
+Tp  = 900.0       # Period of oscillation
+a   = 1.0         # Amplitude at origin
+g   = 9.81        # Gravity
+
+# Bessel functions
+def j0(x):
+    return jn(0.0, x)
+
+def j1(x):
+    return jn(1.0, x)
+
+#DIMENSIONLESS PARAMETERS
+eps = a/h0
+T = Tp*sqrt(g*h0)/L
+A = eps/j0(4.0*pi/T)
+
+# structured mesh
+#points, vertices, boundary = anuga.rectangular_cross(int(1.1*L/dx), int(W/dy), 1.1*L, W, (-1.1*L/2.0, -W/2.0))
+points, vertices, boundary = anuga.rectangular_cross(int(1.1*L/dx), int(W/dy), 1.1*L, W, (0.0, 0.0))
+
+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)
+
+def elevation(x,y):
+    N = len(x)    
+    z = zeros(N, float)
+    for i in range(N):
+        z[i] = (h0/L)*x[i] - h0
+    return z
+domain.set_quantity('elevation', elevation)
+
+def height(x,y):
+    N = len(x)    
+    h = zeros(N, float)
+    for i in range(N):
+        h[i] = h0 - (h0/L)*x[i]
+        if h[i] < 0.0:
+            h[i] = 0.0
+    return h
+domain.set_quantity('height', height)
+
+def stage(x,y):
+    h = height(x,y)
+    z = elevation(x,y)
+    return h+z
+domain.set_quantity('stage', stage)
+
+
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+##def shore(t):
+##    def g(u):
+##        return u + 2.0*A*pi/T*sin(2.0*pi/T*(t+u))    
+##    u = fsolve(g,0.0)
+##    xi = -0.5*u*u + A*cos(2.0*pi/T*(t+u))
+##    position = 1.0 + xi
+##    return position, u           # dimensionless
+
+def prescribe(x,t): 
+    q = zeros(2)
+    def fun(q):                  # Here q=(w, u)
+        f = zeros(2)
+        f[0] = q[0] + 0.5*q[1]**2.0 - A*j0(4.0*pi/T*(1.0+q[0]-x)**0.5)*cos(2.0*pi/T*(t+q[1]))
+        f[1] = q[1] + A*j1(4.0*pi/T*(1.0+q[0]-x)**0.5)*sin(2.0*pi/T*(t+q[1]))/(1+q[0]-x)**0.5
+        return f
+    q = fsolve(fun,q)    
+    return q[0], q[1]            # dimensionless
+
+def f_CG(t): 
+    timing = t*sqrt(g*h0)/L      # dimensionless
+    w, u = prescribe(0.0,timing) # dimensionless
+    wO = w*h0                    # dimensional
+    uO = u*sqrt(g*h0)            # dimensional
+    zO = -h0                     # dimensional
+    hO = wO - zO                 # dimensional
+    pO = uO * hO                 # dimensional
+    #[    'stage', 'Xmomentum', 'Ymomentum']    
+    return [wO,  pO, 0.0]        # dimensional
+
+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
+BTime = anuga.Time_boundary(domain,f_CG)
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': BTime, 'right': Bt, '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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = Tp/48., finaltime = 7000.):
+    #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/carrier_greenspan_periodic/plot_results_cross_section.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/plot_results_cross_section.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/plot_results_cross_section.py	(revision 8734)
@@ -0,0 +1,95 @@
+"""
+    Quick plot of the Carrier-Greenspan outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_carrier_greenspan import *
+
+p_st = util.get_output('carrier_greenspan.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+W1, P1, Z1, H1, U1 = analytic_cg(p2_st.x[v2], p2_st.time[288], h0=5e2, L=5e4, a=1.0, Tp=900.0)
+W2, P2, Z2, H2, U2 = analytic_cg(p2_st.x[v2], p2_st.time[296], h0=5e2, L=5e4, a=1.0, Tp=900.0)
+W3, P3, Z3, H3, U3 = analytic_cg(p2_st.x[v2], p2_st.time[304], h0=5e2, L=5e4, a=1.0, Tp=900.0)
+W4, P4, Z4, H4, U4 = analytic_cg(p2_st.x[v2], p2_st.time[312], h0=5e2, L=5e4, a=1.0, Tp=900.0)
+W5, P5, Z5, H5, U5 = analytic_cg(p2_st.x[v2], p2_st.time[320], h0=5e2, L=5e4, a=1.0, Tp=900.0)
+W6, P6, Z6, H6, U6 = analytic_cg(p2_st.x[v2], p2_st.time[328], h0=5e2, L=5e4, a=1.0, Tp=900.0)
+
+#Plot the stages##############################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[288,v2], 'b.', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], p2_st.stage[296,v2], 'b.')                    # 1*T/6
+pyplot.plot(p2_st.x[v2], p2_st.stage[304,v2], 'b.')                     # 2*T/6
+pyplot.plot(p2_st.x[v2], p2_st.stage[312,v2], 'b.')                     # 3*T/6
+pyplot.plot(p2_st.x[v2], p2_st.stage[320,v2], 'b.')                     # 4*T/6
+pyplot.plot(p2_st.x[v2], p2_st.stage[328,v2], 'b.')  # 5*T/6
+
+pyplot.plot(p2_st.x[v2], W1,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], W2,'r-')
+pyplot.plot(p2_st.x[v2], W3,'r-')
+pyplot.plot(p2_st.x[v2], W4,'r-')
+pyplot.plot(p2_st.x[v2], W5,'r-')
+pyplot.plot(p2_st.x[v2], W6,'r-')
+
+
+pyplot.title('Stage at several instants in time')
+pyplot.ylim(-5.0,5.0)
+pyplot.legend(loc=2)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png')
+
+
+#Plot the momentums##########################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[288,v2], 'b.', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xmom[296,v2], 'b.')                    # 1*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xmom[304,v2], 'b.')                     # 2*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xmom[312,v2], 'b.')                     # 3*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xmom[320,v2], 'b.')                     # 4*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xmom[328,v2], 'b.')  # 5*T/6
+
+pyplot.plot(p2_st.x[v2], P1,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], P2,'r-')
+pyplot.plot(p2_st.x[v2], P3,'r-')
+pyplot.plot(p2_st.x[v2], P4,'r-')
+pyplot.plot(p2_st.x[v2], P5,'r-')
+pyplot.plot(p2_st.x[v2], P6,'r-')
+
+pyplot.title('Xmomentum at several instants in time')
+pyplot.legend(loc=4)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('xmom_plot.png')
+
+
+
+#Plot the velocities#########################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[288,v2], 'b.', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xvel[296,v2], 'b.')                    # 1*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xvel[304,v2], 'b.')                     # 2*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xvel[312,v2], 'b.')                     # 3*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xvel[320,v2], 'b.')                     # 4*T/6
+pyplot.plot(p2_st.x[v2], p2_st.xvel[328,v2], 'b.')  # 5*T/6
+
+pyplot.plot(p2_st.x[v2], U1,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], U2,'r-')
+pyplot.plot(p2_st.x[v2], U3,'r-')
+pyplot.plot(p2_st.x[v2], U4,'r-')
+pyplot.plot(p2_st.x[v2], U5,'r-')
+pyplot.plot(p2_st.x[v2], U6,'r-')
+
+pyplot.title('Xvelocity at several instants in time')
+pyplot.legend(loc=3)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/plot_results_origin_wrt_time.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/plot_results_origin_wrt_time.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/plot_results_origin_wrt_time.py	(revision 8734)
@@ -0,0 +1,23 @@
+from anuga.utilities import plot_utils as util
+from matplotlib import pyplot as pyplot
+
+p2=util.get_output('carrier_greenspan.sww',0.001)
+p=util.get_centroids(p2,velocity_extrapolation=True)
+
+maxx=p.x.max()
+
+minx=p.x.min()
+
+v1 = abs(p.x -minx).argmin()
+#v2 = abs(p.x -0.5*(minx+maxx)).argmin()
+#v3 = abs(p.x -maxx).argmin()
+
+pyplot.plot(p.time, p.stage[:,v1], label='Left edge of domain')
+#pyplot.plot(p.time, p.stage[:,v2], label='Middle of domain')
+#pyplot.plot(p.time, p.stage[:,v3], label='Right edge of domain')
+pyplot.ylim((-1.0,1.0))
+pyplot.legend()
+pyplot.xlabel('Time')
+pyplot.ylabel('Perturbation')
+pyplot.title('Vertical perturbation at the origin over time')
+pyplot.savefig('perturbation_at_origin.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/produce_results.py	(revision 8734)
@@ -0,0 +1,19 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_carrier_greenspan.py')
+    run_validation_script('plot_results_cross_section.py')
+    run_validation_script('plot_results_origin_wrt_time.py')    
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/results.tex	(revision 8734)
@@ -0,0 +1,41 @@
+
+\section{Carrier-Greenspan periodic solution}
+
+The Carrier-Greenspan periodic solution for flows on a sloping beach. References: [1]. Carrier and Greenspan, Journal of Fluid Mechanics, 1958 and [2]. Mungkasi and Roberts, International Journal for Numerical Methods in Fluids, 2012.
+
+\subsection{Results}
+
+
+We should see excellent agreement between the analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{perturbation_at_origin.png}
+\end{center}
+\caption{Velocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/rootsearch.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/rootsearch.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/rootsearch.py	(revision 8734)
@@ -0,0 +1,18 @@
+## module rootsearch
+''' x1,x2 = rootsearch(f,a,b,dx).
+    Searches the interval (a,b) in increments dx for
+    the bounds (x1,x2) of the smallest root of f(x). 
+    Returns x1 = x2 = None if no roots were detected.
+
+    Taken from the book Numerical Methods in Engineering with Python
+    by J. Kiusalaas    
+'''   
+def rootsearch(f,a,b,dx):
+    x1 = a; f1 = f(a)
+    x2 = a + dx; f2 = f(x2)
+    while f1*f2 > 0.0:
+        if x1  >=  b: return None,None
+        x1 = x2; f1 = f2
+        x2 = x1 + dx; f2 = f(x2)
+    else:
+        return x1,x2
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/swap.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/swap.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/swap.py	(revision 8734)
@@ -0,0 +1,30 @@
+## module swap
+''' swapRows(v,i,j).
+    Swaps rows i and j of vector or matrix [v].
+
+    swapCols(v,i,j).
+    Swaps columns i and j of matrix [v].
+
+def swapRows(v,i,j):
+    if len(v.getshape()) == 1: v[i],v[j] = v[j],v[i]
+    else:
+        temp = v[i].copy()
+        v[i] = v[j]
+        v[j] = temp
+
+Taken from the book Numerical Methods in Engineering with Python
+by J. Kiusalaas        
+'''
+
+def swapRows(v,i,j):
+    if len(v.flat) == 2:
+        v[i],v[j] = v[j],v[i]
+    else:
+        temp = v[i].copy()
+        v[i] = v[j]
+        v[j] = temp      
+
+def swapCols(v,i,j):
+    temp = v[:,j].copy()
+    v[:,j] = v[:,i]
+    v[:,i] = temp
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/util.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/util.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_periodic/util.py	(revision 8734)
@@ -0,0 +1,174 @@
+"""
+This module contains various auxiliary function used by pyvolution.
+"""
+
+def mean(x):
+    from numpy import sum
+    return sum(x)/len(x)
+
+
+def gradient(x0, x1, q0, q1):
+
+    if q1-q0 != 0:
+        a = (q1-q0)/(x1-x0)
+    else:
+        a = 0
+        
+    return a
+
+def  minmod(beta_p,beta_m):
+    if (abs(beta_p) < abs(beta_m)) & (beta_p*beta_m > 0.0):
+        phi = beta_p
+    elif (abs(beta_m) < abs(beta_p)) & (beta_p*beta_m > 0.0):
+        phi = beta_m
+    else:
+        phi = 0.0
+    return phi
+
+def  minmod_kurganov(a,b,c):
+    from numpy import sign
+    if sign(a)==sign(b)==sign(c):
+        return sign(a)*min(abs(a),abs(b),abs(c))
+    else:
+        return 0.0
+
+def  maxmod(a,b):
+    if (abs(a) > abs(b)) & (a*b > 0.0):
+        phi = a
+    elif (abs(b) > abs(a)) & (a*b > 0.0):
+        phi = b
+    else:
+        phi = 0.0
+    return phi
+
+def vanleer(a,b):
+    if abs(a)+abs(b) > 1e-12:
+        return (a*abs(b)+abs(a)*b)/(abs(a)+abs(b))
+    else:
+        return 0.0
+
+def vanalbada(a,b):
+    if a*a+b*b > 1e-12:
+        return (a*a*b+a*b*b)/(a*a+b*b)
+    else:
+        return 0.0
+
+def calculate_wetted_area(x1,x2,z1,z2,w1,w2):
+    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
+        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
+        A = 0.5*(w1-z1)*(x-x1)
+        L = x-x1
+    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
+        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
+        A = 0.5*(w2-z2)*(x2-x)
+        L = x2-x
+    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
+        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
+        A = 0.5*(w2-z2)*(x2-x)
+        L = x2-x
+    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
+        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
+        A = 0.5*(w1-z1)*(x-x1)
+        L = x-x1
+    elif (w1 <= z1) & (w2 <= z2):
+        A = 0.0
+    elif (w1 == z1) & (w2 > z2) & (z2 < z1):
+        A = 0.5*(x2-x1)*(w2-z2)
+    elif (w2 == z2) & (w1 > z1) & (z1 < z2):
+        A = 0.5*(x2-x1)*(w1-z1)
+    return A
+
+
+def calculate_new_wet_area(x1,x2,z1,z2,A):
+    from numpy import sqrt
+    min_centroid_height = 1.0e-3
+    # Assumes reconstructed stage flat in a wetted cell
+    M = (z2-z1)/(x2-x1)
+    L = (x2-x1)
+    min_area = min_centroid_height*L
+    max_area = 0.5*(x2-x1)*abs(z2-z1)
+    if A < max_area:
+        if (z1 < z2):
+            x = sqrt(2*A/M)+x1
+            wet_len = x-x1
+            wc = z1 + sqrt(M*2*A)
+        elif (z2 < z1):
+            x = -sqrt(-2*A/M)+x2
+            wet_len = x2-x 
+            wc = z2+sqrt(-M*2*A)
+        else:
+            wc = A/L+0.5*(z1+z2)
+            wet_len = x2-x1
+    else:
+        wc = 0.5*(z1+z2)+A/L
+        wet_len = x2-x1
+            
+    return wc,wet_len
+
+def calculate_new_wet_area_analytic(x1,x2,z1,z2,A,t):
+    min_centroid_height = 1.0e-3
+    # Assumes reconstructed stage flat in a wetted cell
+    M = (z2-z1)/(x2-x1)
+    L = (x2-x1)
+    min_area = min_centroid_height*L
+    max_area = 0.5*(x2-x1)*abs(z2-z1)
+    w1,uh1 = analytic_cannal(x1,t)
+    w2,uh2 = analytic_cannal(x2,t)
+    if (w1 > z1) & (w2 < z2) & (z1 <= z2):
+        print "test1"
+        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
+        wet_len = x-x1
+    elif (w1 < z1) & (w2 > z2) & (z1 < z2):
+        print "test2"
+        x = ((w2-z1)*(x2-x1)+x1*(z2-z1)-x2*(w2-w1))/(z2-z1+w1-w2)
+        wet_len = x2-x
+    elif (w1 < z1) & (w2 > z2) & (z1 >= z2):
+        print "test3"
+        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
+        wet_len = x2-x
+    elif (w1 > z1) & (w2 < z2) & (z1 > z2):
+        print "test4"
+        x = ((w1-z2)*(x2-x1)+x2*(z2-z1)-x1*(w2-w1))/(z2-z1+w1-w2)
+        wet_len = x-x1
+    elif (w1 >= z1) & (w2 >= z2):
+        print "test5"
+        wet_len = x2-x1 
+    else: #(w1 <= z1) & (w2 <= z2)
+        print "test5"
+        if (w1 > z1) | (w2 > z2):
+            print "ERROR"
+        wet_len = x2-x1        
+    return w1,w2,wet_len,uh1,uh2
+
+def analytic_cannal(C,t):
+    from numpy import zeros,sqrt,sin,cos
+
+    
+    #N = len(C)
+    #u = zeros(N)    ## water velocity
+    #h = zeros(N)    ## water depth
+    x = C
+    g = 9.81
+
+
+    ## Define Basin Bathymetry
+    #z_b = zeros(N) ## elevation of basin
+    #z = zeros(N)   ## elevation of water surface
+    z_infty = 10.0       ## max equilibrium water depth at lowest point.
+    L_x = 2500.0         ## width of channel
+
+    A0 = 0.5*L_x                  ## determines amplitudes of oscillations
+    omega = sqrt(2*g*z_infty)/L_x ## angular frequency of osccilation
+
+    x1 = A0*cos(omega*t)-L_x # left shoreline
+    x2 = A0*cos(omega*t)+L_x # right shoreline
+    if (x >=x1) & (x <= x2):
+        z_b = z_infty*(x**2/L_x**2) ## or A0*cos(omega*t)\pmL_x
+        u = -A0*omega*sin(omega*t)
+        z = z_infty+2*A0*z_infty/L_x*cos(omega*t)*(x/L_x-0.5*A0/(L_x)*cos(omega*t))
+    else:
+       z_b = z_infty*(x**2/L_x**2)
+       u=0.0
+       z = z_b
+    h = z-z_b
+    return z,u*h
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/analytical_cg_transient.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/analytical_cg_transient.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/analytical_cg_transient.py	(revision 8734)
@@ -0,0 +1,81 @@
+"""
+Transient water flows using ANUGA,
+where water driven up a linear sloping beach and time varying boundary.
+Ref: Carrier and Greenspan, Journal of Fluid Mechanics, 1958
+"""
+from scipy import linspace, zeros, array, dot
+from math import log, ceil, sqrt
+from scipy.optimize import fsolve
+
+eps = 0.2
+a_const = 1.5*(1.0+0.9*eps)**0.5
+
+def shoreline(t):
+    def f(u):
+        numer = 5.0*(2.0/a_const*(u+t))**3.0 - (2.0/a_const*(u+t))**5.0
+        denom = (1.0 + 4.0/a_const**2.0 * (u+t)**2.0)**4.0
+        return u - 8.0*eps/a_const*numer/denom    
+    u = fsolve(f,0.0)  # t must be involved here.
+    lam = 2.0/a_const*(u+t)
+    x = -0.5*u**2.0 + eps - eps*(1.0 + 3*lam**2.0 - 2.0*lam**4.0)/(1.0+lam**2.0)**3.0
+    return x, u
+
+def bed(x):
+    return x
+
+def analytical_sol(X,t)    :
+    p, v = shoreline(t)
+    N_X = len(X)
+    Vel = zeros(N_X)
+    Stage = zeros(N_X)
+    q = zeros(2)
+    #q[0] = p
+    #q[1] = v
+    z = bed(X)
+
+    def fun(q): #Here q=(eta, u)
+        fun = zeros(2)
+        fun[0] = q[0] + 0.5*q[1]**2.0 - eps*(1.0 - 2.0*(1.25-(2.0/a_const*(q[1]+t))*1.0j) / ((1.0-(2.0/a_const*(q[1]+t))*1.0j)**2.0 + (16.0*(q[0]-x)/a_const**2.0))**1.5 + 1.5*(1.0-(2.0/a_const*(q[1]+t))*1.0j)**2.0/((1.0-(2.0/a_const*(q[1]+t))*1.0j)**2.0+(16.0*(q[0]-x)/a_const**2.0))**2.5).real
+        fun[1] = q[1] - 8.0*eps/a_const*(1.0/((1.0-(2.0/a_const*(q[1]+t))*1.0j)**2.0 + (16.0*(q[0]-x)/a_const**2.0))**1.5 - 0.75*(1.0-(2.0/a_const*(q[1]+t))*1.0j)/((1.0-(2.0/a_const*(q[1]+t))*1.0j)**2.0+(16.0*(q[0]-x)/a_const**2.0))**2.5).imag
+        return fun
+
+    for m in range(N_X):
+        k = N_X - m-1
+        x = X[k]
+        if X[k] > p:
+            Vel[k] = 0.0
+            Stage[k] = z[k]
+        else:
+            q = fsolve(fun,q)
+            Stage[k] = q[0]
+            Vel[k] = q[1]
+    return Stage, z, Vel
+
+"""
+T = linspace(0.0, 40.0, 2)  # Given time
+X = linspace(-0.5, 0.3, 201)
+for i in range(len(T)):
+    t=T[i]
+    w, z, u = analytical_sol(X,t)
+    from pylab import clf,plot,title,xlabel,ylabel,legend,savefig,show,hold,subplot
+    hold(False)
+    clf()
+    plot1 = subplot(211)
+    plot(X,w, X,z)
+    xlabel('Position')
+    ylabel('Stage')
+    #plot1.set_ylim([0.0,0.25])
+    legend(('Stage', ' Bed'),
+           'lower right', shadow=False)        
+    plot2 = subplot(212)
+    plot(X,u)
+    xlabel('Position')
+    ylabel('Velocity')
+    legend(('Velocity', ' '),
+           'upper right', shadow=False)
+    filename = "cg_"
+    filename += str(i)
+    filename += ".png"
+    savefig(filename)
+    show()
+"""
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/numerical_cg_transient.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/numerical_cg_transient.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/numerical_cg_transient.py	(revision 8734)
@@ -0,0 +1,120 @@
+"""
+Transient water flows  using ANUGA,
+where water driven up a linear sloping beach and time varying boundary.
+Ref: Carrier and Greenspan, Journal of Fluid Mechanics, 1958
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+from numpy import zeros, array, float
+from time import localtime, strftime, gmtime
+from scipy.optimize import fsolve
+from math import sin, pi, exp, sqrt, cos
+from analytical_cg_transient 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 = 'carrier_greenspan_'+time
+output_dir = '.'
+output_file = 'carrier_greenspan'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+#DIMENSIONAL PARAMETERS
+dx  = 500.
+dy  = dx
+L   = 5e4         # Length of channel (m)
+W   = 3*dx        # Width of channel (m)       
+h0 = 5e2          # Height at origin when the water is still
+g   = 9.81        # Gravity
+
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(0.8*L/dx), int(W/dy), 0.8*L, W, (-0.5*L, -0.5*W))
+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)
+
+def elevation(x,y):
+    return h0*x/L
+domain.set_quantity('elevation', elevation)
+
+def stage(x,y):
+    w,z,u = analytical_sol(x/L,0.0)
+    return h0*w
+domain.set_quantity('stage', stage)
+
+
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+def f_CG(t): 
+    timing = t*sqrt(g*h0)/L      # dimensionless
+    w, z, u = analytical_sol(array([-0.5]),timing) # dimensionless
+    wB = w*h0                # dimensional
+    uB = u*sqrt(g*h0)            # dimensional
+    zB = z*h0                     # dimensional
+    hB = wB - zB                 # dimensional
+    pB = uB * hB                 # dimensional
+    #[    'stage', 'Xmomentum', 'Ymomentum']    
+    return [wB,  pB, 0.0]        # dimensional
+
+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
+BTime = anuga.Time_boundary(domain,f_CG)
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': BTime, '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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 1000., finaltime = 30000.):
+    #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/carrier_greenspan_transient/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/plot_results.py	(revision 8734)
@@ -0,0 +1,91 @@
+"""
+    Quick plot of the Carrier-Greenspan outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_cg_transient import *
+
+p_st = util.get_output('carrier_greenspan.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+#Dimensional parameters
+L   = 5e4         # Length of channel (m)      
+h0  = 5e2          # Height at origin when the water is still
+g   = 9.81        # Gravity
+
+#Dimensionless solution
+W1, Z1, U1 = analytical_sol(p2_st.x[v2]/L, p2_st.time[0]*sqrt(g*h0)/L)
+W2, Z2, U2 = analytical_sol(p2_st.x[v2]/L, p2_st.time[1]*sqrt(g*h0)/L)
+W3, Z3, U3 = analytical_sol(p2_st.x[v2]/L, p2_st.time[30]*sqrt(g*h0)/L)
+
+#Dimensional solution
+W1 = W1*h0                # dimensional
+U1 = U1*sqrt(g*h0)        # dimensional
+Z1 = Z1*h0                # dimensional
+H1 = W1 - Z1              # dimensional
+P1 = U1 * H1              # dimensional
+
+W2 = W2*h0                # dimensional
+U2 = U2*sqrt(g*h0)        # dimensional
+Z2 = Z2*h0                # dimensional
+H2 = W2 - Z2              # dimensional
+P2 = U2 * H2              # dimensional
+
+W3 = W3*h0                # dimensional
+U3 = U3*sqrt(g*h0)        # dimensional
+Z3 = Z3*h0                # dimensional
+H3 = W3 - Z3              # dimensional
+P3 = U3 * H3              # dimensional
+
+#Plot the stages##############################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[0,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.stage[1,v2], 'b.') 
+pyplot.plot(p2_st.x[v2], p2_st.stage[30,v2], 'b.')                    
+pyplot.plot(p2_st.x[v2], W1,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], W2,'r-')
+pyplot.plot(p2_st.x[v2], W3,'r-')
+pyplot.title('Stage at several instants in time')
+#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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[0,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[1,v2], 'b.') 
+pyplot.plot(p2_st.x[v2], p2_st.xmom[30,v2], 'b.')                    
+pyplot.plot(p2_st.x[v2], P1,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], P2,'r-')
+pyplot.plot(p2_st.x[v2], P3,'r-')
+pyplot.title('Xmomentum at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('xmom_plot.png')
+
+
+#Plot the velocities#########################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[0,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[1,v2], 'b.')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[30,v2], 'b.')                    
+pyplot.plot(p2_st.x[v2], U1,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], U2,'r-')
+pyplot.plot(p2_st.x[v2], U3,'r-')
+pyplot.title('Xvelocity at several instants in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_cg_transient.py')
+    run_validation_script('plot_results.py')   
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/carrier_greenspan_transient/results.tex	(revision 8734)
@@ -0,0 +1,32 @@
+
+\section{Carrier-Greenspan transient solution}
+
+The Carrier-Greenspan transient solution for flows on a sloping beach. References: Carrier and Greenspan, Journal of Fluid Mechanics, 1958.
+
+\subsection{Results}
+
+
+We should see excellent agreement between the analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/analytical_dam_break_dry.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/analytical_dam_break_dry.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/analytical_dam_break_dry.py	(revision 8734)
@@ -0,0 +1,174 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+#       untitled.py
+#       
+#       Copyright 2011 Stephen Roberts <steve@babbage>
+#       
+#       This program is free software; you can redistribute it and/or modify
+#       it under the terms of the GNU General Public License as published by
+#       the Free Software Foundation; either version 2 of the License, or
+#       (at your option) any later version.
+#       
+#       This program is distributed in the hope that it will be useful,
+#       but WITHOUT ANY WARRANTY; without even the implied warranty of
+#       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#       GNU General Public License for more details.
+#       
+#       You should have received a copy of the GNU General Public License
+#       along with this program; if not, write to the Free Software
+#       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+#       MA 02110-1301, USA.
+
+    
+def wu(h2, h0=1.0, h1=10.0):
+    """ The intermediate state in dambreak problem
+    must satisfy wu(h2) = 0.0
+    
+    Taken from formula (3.14) from Zoppou and Roberts SWW notes
+    """
+
+    import math
+
+    h21 = h2/h1
+    h01 = h0/h1
+
+    #print h0,h1,h01
+    
+
+   # Calculate
+    #h21**3 -9.0*(h21**2)*h01 +16.0*(h21**(1.5))*h01 -h21*h01*(h01+8.0) +h01**3
+
+
+
+    sqrth21 = math.sqrt(h21)
+
+    #calc_1 = h21**3 -9.0*(h21**2)*h01 +16.0*(h21**(1.5))*h01 -h21*h01*(h01+8.0) +h01**3
+
+    calc_2 = h21*(sqrth21*(sqrth21*(h21-9*h01)+16.0*h01)- h01*(h01+8.0)) + h01**3
+
+
+    #print calc_1, calc_2
+    
+    return calc_2
+
+
+def calc_h2(h0=1.0, h1=10.0):
+
+    assert h1 >  0.0
+    assert h0 <=  h1
+    assert h0 >= 0.0
+
+
+    h01 = h0/h1
+
+    if h01 >= 1 - 1.0e-4:
+        return 0.5*(h0+h1)
+
+
+    from scipy.optimize import brentq
+    return brentq(wu, h0, h1, args=(h0,h1))
+
+
+
+
+
+def dam_break(x,t,h0=1.0,h1 = 10.0):
+    import math
+    from anuga import g
+
+    h2 = calc_h2(h0,h1)
+    u2 = 2.0*(math.sqrt(g*h1) - math.sqrt(g*h2))
+    #u2 = 2.0*(h1-h2)*g/(math.sqrt(g*h1) + math.sqrt(g*h2))
+
+    try:
+        s = u2*h2/(h2 - h0)
+    except ZeroDivisionError:
+        s = math.sqrt(g*h2)
+
+    c1 = math.sqrt(g*h1)
+    c2 = math.sqrt(g*h2)
+
+    # use numpy.select
+    if x < -t*c1:
+        h = h1
+        u = 0.0
+    elif x < t*(u2 - c2):
+        h = 1.0/g*(2.0/3.0 *c1 - 1.0/3.0*x/t)**2
+        u = 2.0/3.0*(c1 + x/t)
+    elif x < s*t:
+        h = h2
+        u = u2
+    else:
+        h = h0
+        u = 0
+
+    return h,u
+
+
+
+
+
+
+def vec_dam_break(x,t,h0=1.0,h1 = 10.0):
+    import math
+    import numpy
+    from anuga import g
+
+    msg = 'Argument x should be a numpy array'
+    assert isinstance(x,numpy.ndarray), msg
+
+    h2 = calc_h2(h0,h1)
+    u2 = 2.0*(math.sqrt(g*h1) - math.sqrt(g*h2))
+
+    try:
+        s = u2*h2/(h2 - h0)
+    except ZeroDivisionError:
+        s = math.sqrt(g*h2)
+
+    c1 = math.sqrt(g*h1)
+    c2 = math.sqrt(g*h2)
+
+
+    condlist = [ x < -t*c1, x < t*(u2-c2), x < s*t]
+    hchoicelist = [h1, 1.0/g*(2.0/3.0 *c1 - 1.0/3.0*x/t)**2, h2 ]
+    uchoicelist = [0.0, 2.0/3.0*(c1 + x/t) , u2 ]
+
+
+    h = numpy.select(condlist,hchoicelist,default=h0)
+    u = numpy.select(condlist,uchoicelist,default=0.0)
+
+    return h,u
+
+
+
+if __name__ == "__main__":
+    import numpy
+    import pylab
+
+    h1 = 10.0
+    h0 = 1.0
+    
+    n = 100
+    x = -1.0 + numpy.arange(0,2*n+1)/float(n)
+    t = 0.075
+
+    try:
+        h,u = vec_dam_break(2.0,t,h0=h0)
+    except:
+        pass
+    else:
+        raise Exception('should have raised an exception')
+    
+    h,u = vec_dam_break(x,t,h0=h0)
+
+    m = n/4
+    hstar, ustar = dam_break(x[m],t,h0=h0)
+
+    assert hstar == h[m]
+    assert ustar == u[m]
+
+    pylab.plot(x,h)
+    #pylab.plot(x,u)
+
+    pylab.show()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/numerical_dam_break_dry.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/numerical_dam_break_dry.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/numerical_dam_break_dry.py	(revision 8734)
@@ -0,0 +1,117 @@
+"""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 anuga
+from anuga import Domain as Domain
+from numpy import zeros, float
+
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+#time = strftime('%Y%m%d_%H%M%S',localtime())
+
+#output_dir = 'dam_break_'+time
+output_dir = '.'
+output_file = 'dam_break'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1.
+dy = dx
+L = 1000.
+W = 5*dx
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (-L/2.0, -W/2.0))
+
+#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('elevation',0.0)
+domain.set_quantity('friction', 0.0)
+
+h0 = 10.0
+h1 = 0.0
+
+def height(x,y):
+    z = zeros(len(x), float)
+    for i in range(len(x)):
+        if x[i]<=0.0:
+            z[i] = h0
+        else:
+            z[i] = h1
+    return z
+domain.set_quantity('stage', height)
+
+#-----------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+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': Bt, 'right': Bt, '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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# 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.5, finaltime = 50.):
+    #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/dam_break_dry/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/plot_results.py	(revision 8734)
@@ -0,0 +1,72 @@
+"""
+    Quick plot of the dam break outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+import analytical_dam_break_dry as analytic
+
+p_st = util.get_output('dam_break.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+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)
+
+h0 = 1e-9
+h1 = 10.0
+
+h10,u10 = analytic.vec_dam_break(p2_st.x[v2], p2_st.time[10], h0=h0, h1=h1)
+h50,u50 = analytic.vec_dam_break(p2_st.x[v2], p2_st.time[50], h0=h0, h1=h1)
+#h100,u100 = analytic.vec_dam_break(p2_st.x[v2], p2_st.time[100], h0=h0, h1=h1)
+
+#Plot stages
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[10,v2],'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.stage[50,v2], 'b.')
+#pyplot.plot(p2_st.x[v2], p2_st.stage[100,v2], 'b.')
+pyplot.plot(p2_st.x[v2], h10,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], h50,'r-')
+#pyplot.plot(p2_st.x[v2], h100,'r-')
+pyplot.title('Stage at several instants in time')
+pyplot.legend(loc=3)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png')
+
+
+#Plot xmomentum
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[10,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[50,v2], 'b.')
+#pyplot.plot(p2_st.x[v2], p2_st.xmom[100,v2],'b.')
+pyplot.plot(p2_st.x[v2], u10*h10,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u50*h50,'r-')
+#pyplot.plot(p2_st.x[v2], u100*h100,'r-')
+pyplot.title('Xmomentum at several instants in time')
+pyplot.legend(loc=2)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('xmom_plot.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[10,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[50,v2], 'b.')
+#pyplot.plot(p2_st.x[v2], p2_st.xvel[100,v2],'b.')
+pyplot.plot(p2_st.x[v2], u10,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u50,'r-')
+#pyplot.plot(p2_st.x[v2], u100,'r-')
+pyplot.title('Xvelocity at several instants in time')
+pyplot.legend(loc=2)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/produce_results.py	(revision 8734)
@@ -0,0 +1,19 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_dam_break_dry.py')
+    run_validation_script('plot_results.py')
+    run('pdflatex', 'report.tex')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/report.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/report.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/report.tex	(revision 8734)
@@ -0,0 +1,39 @@
+\documentclass[11pt,a4paper]{article}
+
+\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 Dry Dam Break}
+\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/dam_break_dry/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_dry/results.tex	(revision 8734)
@@ -0,0 +1,68 @@
+
+\section{Dam break involving a dry area}
+
+The dam break problem involving a dry area was solved analytically by Ritter~\cite{Ritter1892} as well as Stoker~\cite{Stoker1948, Stoker1957}. The analytical solution exhibits a rarefaction fan as a parabolic curve. As water moves, it involves wetting process over the dry area.
+
+The initial condition is
+\begin{equation} \label{eq:db_dry_init}
+u(x,0)=0, ~~v(x,y)=0, ~~\textrm{and}~~
+h(x,0) = \left\{ \begin{array}{ll}
+h_1 & \textrm{if $x < 0$}\\
+0 & \textrm{if $x > 0$}\\
+\end{array} \right.
+\end{equation}
+where $h_1>0$. The topography is a horizontal flat bed.
+
+
+The analytical solution~\cite{Ritter1892, Stoker1948, Stoker1957} at time $t>0$ is
+\begin{equation}
+h(x) = \left\{ \begin{array}{ll}
+h_1 & \textrm{if $x \leq -t \sqrt{gh_1}$}\\
+h_R=\frac{4}{9g}(\sqrt{gh_1}-\frac{x}{2t})^2 & \textrm{if $-t \sqrt{gh_1} <x \leq 2t\sqrt{gh_1}$}\\
+0 & \textrm{if $x \geq 2t\sqrt{gh_1}$}\\
+\end{array} \right.
+\end{equation}
+which is the free surface and
+\begin{equation}
+u(x) = \left\{ \begin{array}{ll}
+0 & \textrm{if $x \leq -t \sqrt{gh_1}$}\\
+u_R=\frac{2}{3}(\sqrt{gh_1}+\frac{x}{t}) & \textrm{if $-t \sqrt{gh_1} <x \leq 2t\sqrt{gh_1}$}\\
+0 & \textrm{if $x \geq 2t\sqrt{gh_1}$}\\
+\end{array} \right.
+\end{equation}
+which is the velocity.
+
+
+
+\subsection{Results}
+For our test, we consider $h_1=10$ in (\ref{eq:db_dry_init}).
+The following figures show the stage, $x$-momentum, and $x$-velocity at several instants of time. We should see excellent agreement between the analytical and numerical solutions. The wet/dry interface is difficult to resolve and it usually produces large errors.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.8\textwidth]{stage_plot.png}
+%\label{fig:db_dry_stage}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.8\textwidth]{xmom_plot.png}
+%\label{fig:db_dry_xmom}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.8\textwidth]{xvel_plot.png}
+%\label{fig:db_dry_xvel}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/analytical_dam_break_wet.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/analytical_dam_break_wet.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/analytical_dam_break_wet.py	(revision 8734)
@@ -0,0 +1,174 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+#       untitled.py
+#       
+#       Copyright 2011 Stephen Roberts <steve@babbage>
+#       
+#       This program is free software; you can redistribute it and/or modify
+#       it under the terms of the GNU General Public License as published by
+#       the Free Software Foundation; either version 2 of the License, or
+#       (at your option) any later version.
+#       
+#       This program is distributed in the hope that it will be useful,
+#       but WITHOUT ANY WARRANTY; without even the implied warranty of
+#       MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#       GNU General Public License for more details.
+#       
+#       You should have received a copy of the GNU General Public License
+#       along with this program; if not, write to the Free Software
+#       Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+#       MA 02110-1301, USA.
+
+    
+def wu(h2, h0=1.0, h1=10.0):
+    """ The intermediate state in dambreak problem
+    must satisfy wu(h2) = 0.0
+    
+    Taken from formula (3.14) from Zoppou and Roberts SWW notes
+    """
+
+    import math
+
+    h21 = h2/h1
+    h01 = h0/h1
+
+    #print h0,h1,h01
+    
+
+   # Calculate
+    #h21**3 -9.0*(h21**2)*h01 +16.0*(h21**(1.5))*h01 -h21*h01*(h01+8.0) +h01**3
+
+
+
+    sqrth21 = math.sqrt(h21)
+
+    #calc_1 = h21**3 -9.0*(h21**2)*h01 +16.0*(h21**(1.5))*h01 -h21*h01*(h01+8.0) +h01**3
+
+    calc_2 = h21*(sqrth21*(sqrth21*(h21-9*h01)+16.0*h01)- h01*(h01+8.0)) + h01**3
+
+
+    #print calc_1, calc_2
+    
+    return calc_2
+
+
+def calc_h2(h0=1.0, h1=10.0):
+
+    assert h1 >  0.0
+    assert h0 <=  h1
+    assert h0 >= 0.0
+
+
+    h01 = h0/h1
+
+    if h01 >= 1 - 1.0e-4:
+        return 0.5*(h0+h1)
+
+
+    from scipy.optimize import brentq
+    return brentq(wu, h0, h1, args=(h0,h1))
+
+
+
+
+
+def dam_break(x,t,h0=1.0,h1 = 10.0):
+    import math
+    from anuga import g
+
+    h2 = calc_h2(h0,h1)
+    u2 = 2.0*(math.sqrt(g*h1) - math.sqrt(g*h2))
+    #u2 = 2.0*(h1-h2)*g/(math.sqrt(g*h1) + math.sqrt(g*h2))
+
+    try:
+        s = u2*h2/(h2 - h0)
+    except ZeroDivisionError:
+        s = math.sqrt(g*h2)
+
+    c1 = math.sqrt(g*h1)
+    c2 = math.sqrt(g*h2)
+
+    # use numpy.select
+    if x < -t*c1:
+        h = h1
+        u = 0.0
+    elif x < t*(u2 - c2):
+        h = 1.0/g*(2.0/3.0 *c1 - 1.0/3.0*x/t)**2
+        u = 2.0/3.0*(c1 + x/t)
+    elif x < s*t:
+        h = h2
+        u = u2
+    else:
+        h = h0
+        u = 0
+
+    return h,u
+
+
+
+
+
+
+def vec_dam_break(x,t,h0=1.0,h1 = 10.0):
+    import math
+    import numpy
+    from anuga import g
+
+    msg = 'Argument x should be a numpy array'
+    assert isinstance(x,numpy.ndarray), msg
+
+    h2 = calc_h2(h0,h1)
+    u2 = 2.0*(math.sqrt(g*h1) - math.sqrt(g*h2))
+
+    try:
+        s = u2*h2/(h2 - h0)
+    except ZeroDivisionError:
+        s = math.sqrt(g*h2)
+
+    c1 = math.sqrt(g*h1)
+    c2 = math.sqrt(g*h2)
+
+
+    condlist = [ x < -t*c1, x < t*(u2-c2), x < s*t]
+    hchoicelist = [h1, 1.0/g*(2.0/3.0 *c1 - 1.0/3.0*x/t)**2, h2 ]
+    uchoicelist = [0.0, 2.0/3.0*(c1 + x/t) , u2 ]
+
+
+    h = numpy.select(condlist,hchoicelist,default=h0)
+    u = numpy.select(condlist,uchoicelist,default=0.0)
+
+    return h,u
+
+
+
+if __name__ == "__main__":
+    import numpy
+    import pylab
+
+    h1 = 10.0
+    h0 = 1.0
+    
+    n = 100
+    x = -1.0 + numpy.arange(0,2*n+1)/float(n)
+    t = 0.075
+
+    try:
+        h,u = vec_dam_break(2.0,t,h0=h0)
+    except:
+        pass
+    else:
+        raise Exception('should have raised an exception')
+    
+    h,u = vec_dam_break(x,t,h0=h0)
+
+    m = n/4
+    hstar, ustar = dam_break(x[m],t,h0=h0)
+
+    assert hstar == h[m]
+    assert ustar == u[m]
+
+    pylab.plot(x,h)
+    #pylab.plot(x,u)
+
+    pylab.show()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/numerical_dam_break_wet.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/numerical_dam_break_wet.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/numerical_dam_break_wet.py	(revision 8734)
@@ -0,0 +1,122 @@
+"""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, 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 = 'dam_break_'+time
+output_dir = '.'
+output_file = 'dam_break'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1.
+dy = dx
+L = 1000.
+W = 5*dx
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (-L/2.0, -W/2.0))
+
+#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('elevation',0.0)
+domain.set_quantity('friction', 0.0)
+
+h0 = 10.0
+h1 = 1.0
+
+def height(x,y):
+    z = zeros(len(x), float)
+    for i in range(len(x)):
+        if x[i]<=0.0:
+            z[i] = h0
+        else:
+            z[i] = h1
+    return z
+domain.set_quantity('stage', height)
+
+#-----------------------------------------------------------------------------
+# 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': Bt, 'right': Bt, '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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# 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.5, finaltime = 50.):
+    #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/dam_break_wet/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/plot_results.py	(revision 8734)
@@ -0,0 +1,68 @@
+"""
+    Quick plot of the dam break outputs
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+import analytical_dam_break_wet as analytic
+
+p_st = util.get_output('dam_break.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+
+h0 = 1.0
+h1 = 10.0
+
+h10,u10 = analytic.vec_dam_break(p2_st.x[v2], p2_st.time[10], h0=h0, h1=h1)
+h50,u50 = analytic.vec_dam_break(p2_st.x[v2], p2_st.time[50], h0=h0, h1=h1)
+h100,u100 = analytic.vec_dam_break(p2_st.x[v2], p2_st.time[100], h0=h0, h1=h1)
+
+#Plot stages
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[10,v2],'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.stage[50,v2], 'b.')
+pyplot.plot(p2_st.x[v2], p2_st.stage[100,v2], 'b.')
+pyplot.plot(p2_st.x[v2], h10,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], h50,'r-')
+pyplot.plot(p2_st.x[v2], h100,'r-')
+pyplot.title('Stage at several instants in time')
+pyplot.legend(loc=3)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png')
+
+
+#Plot xmomentums
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[10,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[50,v2], 'b.')
+pyplot.plot(p2_st.x[v2], p2_st.xmom[100,v2],'b.')
+pyplot.plot(p2_st.x[v2], u10*h10,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u50*h50,'r-')
+pyplot.plot(p2_st.x[v2], u100*h100,'r-')
+pyplot.title('Xmomentum at several instants in time')
+pyplot.legend(loc=2)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('xmom_plot.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[10,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[50,v2], 'b.')
+pyplot.plot(p2_st.x[v2], p2_st.xvel[100,v2],'b.')
+pyplot.plot(p2_st.x[v2], u10,'r-', label='analytical')
+pyplot.plot(p2_st.x[v2], u50,'r-')
+pyplot.plot(p2_st.x[v2], u100,'r-')
+pyplot.title('Xvelocity at several instants in time')
+pyplot.legend(loc=2)
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/produce_results.py	(revision 8734)
@@ -0,0 +1,19 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_dam_break_wet.py')
+    run_validation_script('plot_results.py')
+    run('pdflatex','report.tex')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/report.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/report.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/report.tex	(revision 8734)
@@ -0,0 +1,38 @@
+\documentclass[11pt,a4paper]{article}
+
+\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 Wet Dam Break}
+\maketitle
+%\tableofcontents
+
+%======================
+%\chapter{Wet 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/dam_break_wet/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/dam_break_wet/results.tex	(revision 8734)
@@ -0,0 +1,70 @@
+
+\section{Dam break on wet areas}
+
+The dam break problem on wet areas was solved analytically by Stoker~\cite{Stoker1948, Stoker1957}. The analytical solution exhibits a rarefaction and involves a shock. Generally this problem is easier to solve numerically than the dry dam break (the dam break on a dry area).
+
+The initial condition is
+\begin{equation} \label{eq:dbp_init_wet}
+u(x,0)=0, ~~v(x,y)=0, ~~\textrm{and}~~
+h(x,0) = \left\{ \begin{array}{ll}
+h_1 & \textrm{if $x < 0$}\\
+h_0 & \textrm{if $x > 0$}\\
+\end{array} \right.
+\end{equation}
+where $h_1>h_0>0$. The topography is a horizontal flat bed.
+
+The analytical solution~\cite{Stoker1948, Stoker1957} to the wet dam break problem is
+\begin{equation}
+h(x) = \left\{ \begin{array}{ll}
+h_1 & \textrm{if $x \leq -t \sqrt{gh_1}$}\\
+h_3=\frac{4}{9g}(\sqrt{gh_1}-\frac{x}{2t})^2 & \textrm{if $-t \sqrt{gh_1} <x \leq t(u_2-\sqrt{gh_2}$})\\
+h_2=\frac{h_0}{2}\bigg(\sqrt{1+\frac{8\dot{\xi}^2}{gh_0}}-1\bigg) & \textrm{if $ t(u_2-\sqrt{gh_2}) <x < t\dot{\xi}$}\\
+h_0 & \textrm{if $x \geq t\dot{\xi}$}\\
+\end{array} \right.
+\end{equation}
+and
+\begin{equation}
+u(x) = \left\{ \begin{array}{ll}
+0 & \textrm{if $x \leq -t \sqrt{gh_1}$}\\
+u_3=\frac{2}{3}(\sqrt{gh_1}+\frac{x}{t}) & \textrm{if $-t \sqrt{gh_1} <x \leq t(u_2-\sqrt{gh_2}$})\\
+u_2=\dot{\xi}-\frac{gh_0}{4\dot{\xi}}\bigg(1+\sqrt{1+\frac{8\dot{\xi}^2}{gh_0}} \bigg) & \textrm{if $ t(u_2-\sqrt{gh_2}) <x < t\dot{\xi}$}\\
+0 & \textrm{if $x \geq t\dot{\xi}$}\\
+\end{array} \right.
+\end{equation}
+at any time $t>0$, where $\dot{\xi}$ is the shock speed constant given by 
+\begin{equation} \label{eq:shock}
+\dot{\xi}=2\sqrt{gh_1}+\frac{gh_0}{4\dot{\xi}}\bigg( 1+\sqrt{1+\frac{8\dot{\xi}^2}{gh_0}}\bigg)-\bigg( 2gh_0 \sqrt{1+\frac{8\dot{\xi}^2}{gh_0}}-2gh_0\bigg)^\frac{1}{2}.
+\end{equation}
+
+
+\subsection{Results}
+
+For our test, we consider $h_1=10$ and $h_0=1$ in (\ref{eq:dbp_init_wet}).
+The following figures show the stage, $x$-momentum, and $x$-velocity at several instants of time.
+We should see excellent agreement between the analytical and numerical solutions. 
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/plot_results.py	(revision 8734)
@@ -0,0 +1,49 @@
+from anuga.utilities import plot_utils as util
+from matplotlib import pyplot as pyplot
+
+p2=util.get_output('data_wave.sww',0.001)
+p=util.get_centroids(p2,velocity_extrapolation=True)
+
+maxx=p.x.max()
+
+minx=p.x.min()
+
+v1 = abs(p.x -minx).argmin()
+v2 = abs(p.x -0.5*(minx+maxx)).argmin()
+v3 = abs(p.x -maxx).argmin()
+
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p.time, p.stage[:,v1], label='Left edge of domain')
+pyplot.plot(p.time, p.stage[:,v2], label='Middle of domain')
+pyplot.plot(p.time, p.stage[:,v3], label='Right edge of domain')
+pyplot.ylim((-1.5,2.0))
+pyplot.legend()
+pyplot.xlabel('Time')
+pyplot.ylabel('Amplitude')
+pyplot.title('Stage over time at 3 points in space')
+pyplot.savefig('wave_atten.png')
+
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p.time, p.xmom[:,v1], label='Left edge of domain')
+pyplot.plot(p.time, p.xmom[:,v2], label='Middle of domain')
+pyplot.plot(p.time, p.xmom[:,v3], label='Right edge of domain')
+#pyplot.ylim((-1.,1.0))
+pyplot.legend()
+pyplot.xlabel('Time')
+pyplot.ylabel('Xmomentum')
+pyplot.title('Xmomentum over time at 3 points in space')
+pyplot.savefig('xmom.png')
+
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p.time, p.ymom[:,v1], label='Left edge of domain')
+pyplot.plot(p.time, p.ymom[:,v2], label='Middle of domain')
+pyplot.plot(p.time, p.ymom[:,v3], label='Right edge of domain')
+#pyplot.ylim((-1.,1.0))
+pyplot.legend()
+pyplot.xlabel('Time')
+pyplot.ylabel('Ymomentum')
+pyplot.title('Ymomentum over time at 3 points in space')
+pyplot.savefig('ymom.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/produce_results.py	(revision 8734)
@@ -0,0 +1,17 @@
+from fabricate import *
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('run_wave.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/results.tex	(revision 8734)
@@ -0,0 +1,35 @@
+\section{Deep Water Wave Propagtion}
+This simulates the free propagation of a sinusoidal wave in deep water. Analytically, the wave should travel through the domain without deformation. However, in early versions of ANUGA, some config settings would lead to the wave dampening out too strongly. This sort of behaviour can have practical implications e.g. for tsunami propagation problems.  
+
+This example can also illustrate difficulties with boundary conditions, as it is difficult to accurately treat the boundary where the wave exits the domain (of course, for this problem, we could do that by exploiting the analytical solution - but this is not possible for general wave propagation problems).
+
+\subsection{Results}
+Figure~\ref{fig:stagewave} shows the time-evolution of the water elevation at three points in the domain. These time series should show the wave propagating without deformation or attenuation (i.e. the wave has the same shape, amplitude, period, mean water level etc at each point).  
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{wave_atten.png}
+\caption{Stage over time at 3 points in space}
+\label{fig:stagewave}
+\end{center}
+\end{figure}
+
+
+The corresponding momentums of Figure~\ref{fig:stagewave} are shown in Figures~\ref{fig:xmom} and~\ref{fig:ymom}.
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom.png}
+\caption{Xmomentum over time at 3 points in space}
+\label{fig:xmom}
+\end{center}
+\end{figure}
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{ymom.png}
+\caption{Ymomentum over time at 3 points in space}
+\label{fig:ymom}
+\end{center}
+\end{figure}
+
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/run_wave.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/run_wave.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/deep_wave/run_wave.py	(revision 8734)
@@ -0,0 +1,135 @@
+"""Simple water flow example using ANUGA
+
+Will Powers example of a simple sinusoidal wave. 
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+import sys
+import anuga
+from anuga import Domain
+
+#from balanced_dev import *
+#from balanced_dev import Domain as Domain
+
+from math import cos
+import numpy as num
+from time import localtime, strftime, gmtime
+from os import sep
+
+
+#--------------------------------
+# Get Default values for basic
+# algorithm parameters.
+#--------------------------------
+from anuga.utilities.argparsing import parse_standard_args
+alg, cfl = parse_standard_args()
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+output_dir = '.'
+output_file = 'data_wave'
+
+#copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+sep)
+
+interactive_visualisation = False
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 500.
+dy = dx
+L = 100000.
+W = 10*dx
+
+# structured mesh
+points, vertices, boundary = anuga.rectangular_cross(int(L/dx), int(W/dy), L, W, (0.0, -W/2))
+
+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('elevation',-100.0)
+domain.set_quantity('friction', 0.00)
+domain.set_quantity('stage', 0.0)            
+
+#-----------------------------------------------------------------------------
+# 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
+#Bz = swb2_boundary_conditions.zero_mass_flux_zero_momentum_boundary(domain) # Strong reflections
+#Bz1 = swb2_boundary_conditions.zero_mass_flux_zero_t_transmissive_n_momentum_boundary(domain) # Strong reflections
+#Bz2 = swb2_boundary_conditions.zero_mass_flux_zero_n_transmissive_t_momentum_boundary(domain) # Strong reflections
+#Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain) # Strong reflections
+#Bd = anuga.Dirichlet_boundary([1,0.,0.]) # Constant boundary values
+amplitude = 1
+wave_length = 300.0
+#Bw = anuga.Time_boundary(domain=domain, function=lambda t: [(-amplitude*sin((1./wave_length)*t*2*pi)), 0.0, 0.0])
+## Sawtooth?
+#                   function=lambda t: [(-8.0*(sin((1./180.)*t*2*pi))+(1./2.)*sin((2./180.)*t*2*pi)+(1./3.)*sin((3./180.)*t*2*pi)), 0.0, 0.0])
+## Sharp rise, linear fall
+#                   function=lambda t: [(5.0*(-((t-0.)/300.)*(t<300.)-cos((t-300.)*2.*pi*(1./240.))*(t>=300. and t<420.)+(1.-(t-420.)/300.)*(t>=420. and t <720.))), 0.0, 0.0])
+#                   function=lambda t: [amplitude*(1.-2.*(pi*(1./720.)*(t-720.))**2)/exp((pi*(1./720.)*(t-720.))**2) , 0.0, 0.0])
+#                   function=lambda t: [(-8.0*sin((1./720.)*t*2*pi))*((t<720.)-0.5*(t<360.)), 0.0, 0.0])
+def waveform(t):
+    return [amplitude*sin((1./wave_length)*t*2*pi),0.0,0.0]
+
+Bw2 = anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, waveform)
+#Bw3 = swb2_boundary_conditions.Transmissive_momentum_nudge_stage_boundary(domain, waveform)
+Bw3 = anuga.Time_boundary(domain, waveform)
+
+def zero_fun(t):
+    return 0.0                    
+#Bw0 = anuga.shallow_water.boundaries.Transmissive_n_momentum_zero_t_momentum_set_stage_boundary(domain, zero_fun)
+#Bw4 = anuga.Transmissive_momentum_set_stage_boundary(domain, waveform)
+
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': Bw2, 'right': Bt, 'top': Br, 'bottom': Br})
+
+
+#===============================================================================
+if interactive_visualisation:
+    from anuga.visualiser import RealtimeVisualiser
+    vis = RealtimeVisualiser(domain)
+    vis.render_quantity_height("stage", zScale =10000, dynamic=True)
+    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+    vis.start()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+
+for t in domain.evolve(yieldstep = 1e1, finaltime = 2e4):
+    domain.write_time()
+    if interactive_visualisation:
+        vis.update()
+
+if interactive_visualisation:
+    vis.evolveFinished()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/numerical_immersed_bump.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/numerical_immersed_bump.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/numerical_immersed_bump.py	(revision 8734)
@@ -0,0 +1,108 @@
+"""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, 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 = 'immersed_bump'+time
+output_dir = '.'
+output_file = 'immersed_bump'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1.
+dy = dx
+L = 25.
+W = 5*dx
+
+# 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) 
+
+#------------------------------------------------------------------------------
+# 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', 0.5)
+
+def bed_elevation(x,y):
+    z = zeros(len(x), float)
+    for i in range(len(x)):
+        if 8.0 < x[i] < 12.0:
+            z[i] = 0.2 - 0.05*(x[i]-10.)**2
+        else:
+            z[i] = 0.0
+    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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 0.1, finaltime = 5.):
+    #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/lake_at_rest_immersed_bump/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/plot_results.py	(revision 8734)
@@ -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
+
+p_st = util.get_output('immersed_bump.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[10,v2],'b.', label='numerical stage')
+pyplot.plot(p2_st.x[v2], p2_st.elev[v2],'k-', label='discretised bed')
+pyplot.plot(p2_st.x[v2], 0.5*ones(len(p2_st.x[v2])),'r-', label='analytical stage')
+pyplot.title('Stage at an instant in time')
+pyplot.legend(loc='best')
+pyplot.ylim([0.0,0.6])
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png')
+
+
+#Plot xmomentum
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[10,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.savefig('xmom_plot.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[10,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.savefig('xvel_plot.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_immersed_bump.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_immersed_bump/results.tex	(revision 8734)
@@ -0,0 +1,35 @@
+
+\section{Lake at rest with an immersed bump}
+
+This is a simple test if the method is well-balanced.
+
+\subsection{Results}
+
+
+We should see excellent agreement between the analytical and numerical solutions if the method is well-balanced.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/numerical_steep_island.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/numerical_steep_island.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/numerical_steep_island.py	(revision 8734)
@@ -0,0 +1,136 @@
+"""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, 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 = 'steep_island'+time
+output_dir = '.'
+output_file = 'steep_island'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 1.
+dy = dx
+L = 2000.
+W = 5*dx
+
+# 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) 
+
+#------------------------------------------------------------------------------
+# 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)
+
+def stage_flat(x,y):
+    w=zeros(len(x))
+    for i in range(len(x)):
+        w[i]=4.5
+    return w
+domain.set_quantity('stage', stage_flat)
+
+def bed_elevation(x,y):
+    z=zeros(len(x))
+    for i in range(len(x)):
+        if 0 <= x[i] < 200.0:
+            z[i] = -0.01*(x[i]-200) + 4.0
+        elif 200.0 <= x[i] < 300.0:
+            z[i] = -0.02*(x[i]-200) + 4.0
+        elif 300.0 <= x[i] < 400.0:
+            z[i] = -0.01*(x[i]-300) + 2.0
+        elif 400.0 <= x[i] < 550.0:
+            z[i] = (-1/75.0)*(x[i]-400.0) + 2.0
+        elif 550.0 <= x[i] < 700.0:
+            z[i] = (1/11250)*(x[i]-550)*(x[i]-550)
+        elif 700.0 <= x[i] < 800.0:
+            z[i] = 0.03*(x[i]-700)
+        elif 800.0 <= x[i] < 900.0:
+            z[i] = -0.03*(x[i]-800) + 3.0
+        elif 900.0 <= x[i] < 1000.0:
+            z[i] = 6.0
+        elif 1000.0 <= x[i] < 1400.0:
+            z[i] = (-1.0/20000)*(x[i]-1000)*(x[i]-1400)
+        elif 1400.0 <= x[i] < 1500.0:
+            z[i] = 0.0
+        elif 1500.0 <= x[i] < 1700.0:
+            z[i] = 3.0
+        elif 1700.0 <= x[i] < 1800.0:
+            z[i] = -0.03*(x[i]-1700) + 3.0
+        else:
+            z[i] = (4.5/40000)*(x[i]-1800)*(x[i]-1800) + 2.0
+    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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 0.1, finaltime = 5.):
+    #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/lake_at_rest_steep_island/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/plot_results.py	(revision 8734)
@@ -0,0 +1,64 @@
+"""
+    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
+
+p_st = util.get_output('steep_island.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+z=zeros(len(p2_st.x[v2]))
+xx = p2_st.x[v2]
+for i in range(len(p2_st.x[v2])):
+    if 0 <= xx[i] <= 200.0:
+        z[i] = max(-0.01*(xx[i]-200) + 4.0, 4.5)
+    elif 900.0 <= xx[i] <= 1000.0:
+        z[i] = 6.0
+    elif 1800.0 <= xx[i] <= 2000.0:
+        z[i] = max((4.5/40000)*(xx[i]-1800)**2 + 2.0, 4.5)
+    else:
+        z[i] = 4.5
+
+#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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[40,v2],'b.', label='numerical stage')
+pyplot.plot(p2_st.x[v2], p2_st.elev[v2],'k-', label='discretised bed')
+pyplot.plot(p2_st.x[v2], z,'r-', label='analytical stage')
+pyplot.title('Stage at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('stage_plot.png')
+
+
+#Plot xmomentum
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[40,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.savefig('xmom_plot.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[40,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.savefig('xvel_plot.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities 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')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/lake_at_rest_steep_island/results.tex	(revision 8734)
@@ -0,0 +1,35 @@
+
+\section{Lake at rest with a steep island}
+
+This is a test if the method is well-balanced. Furthermore, we test if the wet/dry interface has been correctly treated.
+
+\subsection{Results}
+
+
+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}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/analytical_parabolic_basin.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/analytical_parabolic_basin.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/analytical_parabolic_basin.py	(revision 8734)
@@ -0,0 +1,23 @@
+"""
+Analytical solution of Thacker for oscillation on
+a parabolic basin.
+"""
+from numpy import zeros,sqrt,sin,cos
+
+def analytic_cannal(x,t, D0=4., L=10., A=2., g=9.81):
+    omega = sqrt(2*D0*g)/L
+    N = len(x)
+    u = zeros(N)              ## water velocity
+    h = zeros(N)              ## water depth
+    ## Define Basin Bathymetry
+    z = zeros(N)            ## elevation of basin
+    w = zeros(N)              ## elevation of water surface
+    for i in range(N):
+        z[i] = D0*(x[i]**2/L**2)
+        u[i] = -A*omega*sin(omega*t)
+        w[i] = D0 + 2*A*D0/L**2 * cos(omega*t)*( x[i] - A/2.*cos(omega*t))
+        if w[i] <= z[i] :
+            u[i] = 0.0
+            w[i] = z[i]
+    h = w - z
+    return u,h,w,z
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/numerical_parabolic_basin.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/numerical_parabolic_basin.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/numerical_parabolic_basin.py	(revision 8734)
@@ -0,0 +1,108 @@
+"""Parabolic channel oscilation example
+"""
+#---------
+#Import Modules
+#--------
+import anuga
+import numpy
+from math import sqrt, cos, sin, pi
+from numpy import asarray
+from time import localtime, strftime, gmtime
+from anuga.geometry.polygon import inside_polygon, is_inside_triangle
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+timer = strftime('%Y%m%d_%H%M%S',localtime())
+output_dir = '.'  #'parabola_'+timer
+output_file = 'parabola'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+#-------------------------------------------------------------------------------
+# Domain
+#-------------------------------------------------------------------------------
+m = 200
+n = 10
+lenx = 40.0
+leny = 2.0
+origin = (-20.0, -1.0)
+
+points, elements, boundary = anuga.rectangular_cross(m, n, lenx, leny, origin)
+domain = anuga.Domain(points, elements, boundary)
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)
+#domain.set_minimum_allowed_height(0.01)
+
+
+#------------------------------------------------------------------------------
+# 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)
+
+#------------------
+# Define topography
+#------------------
+
+# Parameters for analytical solution
+D0=4.0
+L=10.0
+A = 2.0
+
+def topography(x,y):
+        return  (D0/(L**2.))*x**2.
+
+def stage_init(x,y):
+    wat= D0 + (2.0*A*D0/L**2.)*(x-A/2.0) # Water elevation inside the parabola
+    top=topography(x,y) # Bed elevation
+    # Return the maximum of the water elevation and the bed elvation
+    return wat*(wat>top) + top*(wat<=top)
+
+
+domain.set_quantity('elevation',topography)     # Use function for elevation
+domain.set_quantity('friction',0.0)            # No friction
+domain.set_quantity('stage', stage_init)        # Constant negative initial stage
+
+
+#--------------------------
+# Setup boundary conditions
+#--------------------------
+Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
+#Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
+#Bd=anuga.Dirichlet_boundary([-0.2,0.,0.])       # Constant boundary values -- not used in this example
+#Bw=anuga.Time_boundary(domain=domain,
+#	f=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,0.0,0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
+
+#----------------------------------------------
+# Associate boundary tags with boundary objects
+#----------------------------------------------
+domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom':Br})
+
+
+#---------------------------------------------
+# 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
+
+
+#------------------------------
+#Evolve the system through time
+#------------------------------
+for t in domain.evolve(yieldstep=0.05,finaltime=10.0):
+    print domain.timestepping_statistics()
+print 'Finished'
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/plot_results_cross_section.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/plot_results_cross_section.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/plot_results_cross_section.py	(revision 8734)
@@ -0,0 +1,60 @@
+"""
+    Quick plot of the dam break outputs
+
+"""
+from numpy import zeros
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from math import sqrt, pi, cos, sin
+from analytical_parabolic_basin import analytic_cannal
+
+# Get the sww file
+p_st = util.get_output('parabola.sww')
+p2_st=util.get_centroids(p_st) #(p_st, velocity_extrapolation=True)
+
+index_origin = 3978
+v = p2_st.y[index_origin]
+v2=(p2_st.y==v)
+
+
+# Compute the analytical solution
+time_level = 200
+u,h,w,z = analytic_cannal(p2_st.x[v2], p2_st.time[time_level])
+
+#Plot stages
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[time_level,v2],'b.', label='numerical stage')
+pyplot.plot(p2_st.x[v2], p2_st.elev[v2],'k-', label='bed elevation')
+pyplot.plot(p2_st.x[v2], w,'r-', label='analytical stage')
+pyplot.title('Stage at an instant in time')
+pyplot.legend(loc='upper center')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('cross_section_stage.png')
+
+
+#Plot xmomentum
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[time_level,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], u*h,'r-', label='analytical')
+pyplot.title('Xmomentum at an instant in time')
+pyplot.legend(loc='lower right')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('cross_section_xmom.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[time_level,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], u,'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='upper center')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('cross_section_xvel.png')
+
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/plot_results_origin_wrt_time.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/plot_results_origin_wrt_time.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/plot_results_origin_wrt_time.py	(revision 8734)
@@ -0,0 +1,65 @@
+"""View results of runup.py
+"""
+#---------------
+# Import Modules
+#---------------
+#import anuga
+#import struct
+import numpy
+from matplotlib import pyplot as pyplot
+from anuga.utilities import plot_utils as util
+#from analytical_parabolic_basin import analytic_cannal
+
+p=util.get_output('parabola.sww', 0.01)
+p2=util.get_centroids(p, velocity_extrapolation=True)
+
+### Define some 'lines' along which to plot
+##v=(p2.y==p2.y[0])
+
+# Compute the analytical solution
+D0=4.0
+L=10.0
+A=2.0
+g=9.81
+
+omega=numpy.sqrt(2*D0*g)/L
+T= 2*numpy.pi/omega
+index_origin = 3978
+
+# Free surface over time at the origin of spatial domain.
+w = D0 + 2*A*D0/(L**2)*numpy.cos(omega*p2.time)*( p2.x[index_origin] -A/2.*numpy.cos(omega*p2.time))
+#w = w*(w>p2.elev[index_origin])+p2.elev[index_origin]*(w<=p2.elev[index_origin]) #not needed
+
+# Velocity at the origin of the spatial domain.
+u = -A*omega*numpy.sin(omega*p2.time)
+
+#Plot stages
+pyplot.clf()
+pyplot.plot(p2.time,w, color='blue', label='analytical')
+pyplot.plot(p2.time,p2.stage[:,index_origin], color='green', label='numerical')
+pyplot.legend(loc='best')
+pyplot.title('Stage at the centre of the parabolic basin over time')
+pyplot.xlabel('Time')
+pyplot.ylabel('Stage')
+pyplot.savefig('Stage_centre.png')
+
+#Plot xmomentums
+pyplot.clf()
+pyplot.plot(p2.time,u*w, color='blue', label='analytical') # At the origin of space, w=h as always wet.
+pyplot.plot(p2.time,p2.xmom[:,index_origin], color='green', label='numerical')
+pyplot.legend(loc='best')
+pyplot.title('Xmomentum at the centre of the parabolic basin over time')
+pyplot.xlabel('Time')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('Xmom_centre.png')
+
+#Plot velocities
+pyplot.clf()
+pyplot.plot(p2.time,u, color='blue', label='analytical')
+pyplot.plot(p2.time,p2.xvel[:,index_origin], color='green', label='numerical')
+pyplot.legend(loc='best')
+pyplot.title('Xvelocity at the centre of the parabolic basin over time')
+pyplot.xlabel('Time')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('Xvel_centre.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_parabolic_basin.py')
+    run_validation_script('plot_results_cross_section.py')
+    run_validation_script('plot_results_origin_wrt_time.py')    
+
+def clean():
+    autoclean()
+
+main()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/parabolic_basin/results.tex	(revision 8734)
@@ -0,0 +1,60 @@
+\section{1D Parabolic Basin}
+This simulates flow oscillations in a 1D parabolic basin. The analytical solution is described by Thacker (REF), and is periodic. At any instant in time, the free surface elevation is planar, and the velocity is constant (in wet regions). The scenario includes regular wetting and drying, as the flow oscillates back and forth in the basin. As well as testing the ability of the code to do wetting and drying, it will highlight any numerical energy loss or gain, manifest as an increase or decrease in the magnitude of the flow oscillations over long time periods, compared with the analytical solution. 
+
+\subsection{Results}
+
+Figures~\ref{fig:cross_section_stage}--\ref{fig:cross_section_xvel} show a good agreement between analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{cross_section_stage.png}
+\caption{Stage on a cross section of the basin at time $t=$ ?}
+\label{fig:cross_section_stage}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{cross_section_xmom.png}
+\caption{Xmomentum on a cross section of the basin at time $t=$ ?}
+\label{fig:cross_section_xmom}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{cross_section_xvel.png}
+\caption{Xvelocity on a cross section of the basin at time $t=$ ?}
+\label{fig:cross_section_xvel}
+\end{center}
+\end{figure}
+
+
+There should be a good agreement with the analytical solution, although as time goes on, some deviations may appear.  See Figures~\ref{fig:Stage_centre}--\ref{fig:Xvel_centre}.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{Stage_centre.png}
+\caption{Stage over time in the centre of the parabolic basin.}
+\label{fig:Stage_centre}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{Xmom_centre.png}
+\caption{Xmomentum over time in the centre of the parabolic basin.}
+\label{fig:Xmom_centre}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{Xvel_centre.png}
+\caption{Xvelocity over time in the centre of the parabolic basin.}
+\label{fig:Xvel_centre}
+\end{center}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/analytical_paraboloid_basin.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/analytical_paraboloid_basin.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/analytical_paraboloid_basin.py	(revision 8734)
@@ -0,0 +1,39 @@
+"""
+Thacker's analytical solution for a paraboloid oscillation
+of water on a paraboloid basin
+
+Sudi Mungkasi, ANU 2012
+"""
+from numpy import zeros, sin, cos, sqrt
+
+
+# Set bed elevation
+def bed_elevation(x,y,D0=1000.,L=2500.,R0=2000.,g=9.81):
+    n = x.shape[0]
+    z = 0*x
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        z[i] = -D0*(1.-r*r/L/L)
+    return z
+
+def analytic_sol(x,y,t,D0=1000.,L=2500.,R0=2000.,g=9.81):
+    A = (L**4 - R0**4)/(L**4 + R0**4)
+    omega = 2./L*sqrt(2.*g*D0)
+    z = bed_elevation(x,y)
+    n = x.shape[0]
+    w = 0*x
+    u = 0*x
+    h = 0*x    
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        w[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
+                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
+        u[i] = 0.5*omega*r*A*sin(omega*t) / (1.0-A*cos(omega*t))
+        if x[i] < 0.0:
+            u[i] = -u[i]
+        h[i] = w[i] - z[i]
+        if w[i] < z[i]:
+            w[i] = z[i]
+            u[i] = 0.0
+            h[i] = 0.0
+    return w,u,h
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/numerical_paraboloid_basin.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/numerical_paraboloid_basin.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/numerical_paraboloid_basin.py	(revision 8734)
@@ -0,0 +1,203 @@
+"""Example of shallow water wave equation analytical solution
+consists of a parabolic profile in a parabolic basin. Analytical
+solution to this problem was derived by Thacker
+and used by Yoon and Chou.
+
+   Copyright 2005
+   Christopher Zoppou, Stephen Roberts, ANU, Geoscience Australia
+   Modified by Sudi Mungkasi, ANU, 2011
+"""
+
+#-------------------------------------------------------------------------------
+# Module imports
+#-------------------------------------------------------------------------------
+import sys
+import anuga
+from math import sqrt, cos, sin, pi
+from numpy import asarray
+from time import localtime, strftime, gmtime
+from anuga.geometry.polygon import inside_polygon, is_inside_triangle
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+#timer = strftime('%Y%m%d_%H%M%S',localtime())
+output_dir = '.'  #'paraboloid_'+timer
+output_file = 'paraboloid'
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#-------------------------------------------------------------------------------
+# Domain
+#-------------------------------------------------------------------------------
+n = 75#50
+m = 75#50
+lenx = 8000.0
+leny = 8000.0
+origin = (-4000.0, -4000.0)
+
+points, elements, boundary = anuga.rectangular_cross(m, n, lenx, leny, origin)
+domain = anuga.Domain(points, elements, boundary)
+
+
+#-------------------------------------------------------------------------------
+# Provide file name for storing output
+#-------------------------------------------------------------------------------
+domain.store = True #False
+domain.format = 'sww'
+domain.set_name(output_file)                
+domain.set_datadir(output_dir)
+
+
+#-------------------------------------------------------------------------------
+# Order of scheme
+# Good compromise between
+# limiting and CFL
+#-------------------------------------------------------------------------------
+domain.set_default_order(2)
+domain.set_timestepping_method(2)
+domain.set_beta(1.5)
+domain.set_CFL(0.6)
+domain.smooth = True
+
+
+#-------------------------------------------------------------------------------
+# Decide which quantities are to be stored at each timestep
+#-------------------------------------------------------------------------------
+#domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+
+
+#-------------------------------------------------------------------------------
+# Reduction operation for get_vertex_values
+from anuga.utilities.numerical_tools import mean
+domain.reduction = mean #domain.reduction = min  #Looks better near steep slopes
+
+
+#-------------------------------------------------------------------------------
+# Initial condition
+#-------------------------------------------------------------------------------
+t = 0.0
+D0 = 1000.
+L = 2500.
+R0 = 2000.
+g = 9.81
+
+A = (L**4 - R0**4)/(L**4 + R0**4)
+omega = 2./L*sqrt(2.*g*D0)
+T = pi/omega
+
+#------------------
+# Set bed elevation
+def bed_elevation(x,y):
+    n = x.shape[0]
+    z = 0*x
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        z[i] = -D0*(1.-r*r/L/L)
+    return z
+domain.set_quantity('elevation', bed_elevation)
+
+#----------------------------
+# Set the initial water level
+def stage_init(x,y):
+    z = bed_elevation(x,y)
+    n = x.shape[0]
+    w = 0*x
+    for i in range(n):
+        r = sqrt(x[i]*x[i] + y[i]*y[i])
+        w[i] = D0*((sqrt(1-A*A))/(1.-A*cos(omega*t))
+                -1.-r*r/L/L*((1.-A*A)/((1.-A*cos(omega*t))**2)-1.))
+    if w[i] < z[i]:
+        w[i] = z[i]
+    return w
+domain.set_quantity('stage', stage_init)
+
+#---------
+# Boundary
+#R = anuga.Reflective_boundary(domain)
+#T = anuga.Transmissive_boundary(domain)
+D = anuga.Dirichlet_boundary([0.0, 0.0, 0.0])
+domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D})
+
+#---------------------------------------------
+# 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
+    
+#t = domain.triangles[n_point]
+#print 't = ', t
+#tri = domain.get_vertex_coordinates(n)
+
+filename=domain.get_name()
+file = open(filename,'w')
+
+#----------
+# Evolution
+import time
+t0 = time.time()
+
+time_array = []
+stage_array = []
+Stage     = domain.quantities['stage']
+Xmomentum = domain.quantities['xmomentum']
+Ymomentum = domain.quantities['ymomentum']
+
+#interactive_visualisation = True
+#===============================================================================
+##if interactive_visualisation:
+##    from anuga.visualiser import RealtimeVisualiser
+##    vis = RealtimeVisualiser(domain)
+##    vis.render_quantity_height("stage", zScale =1000, dynamic=True)
+##    vis.colour_height_quantity('stage', (1.0, 0.5, 0.5))
+##    vis.start()
+#===============================================================================
+
+for t in domain.evolve(yieldstep = 1.0, finaltime = 200.0):
+    domain.write_time()
+
+    #tri_array = asarray(tri)
+    #t_array = asarray([[0,1,2]])
+    #interp = Interpolation(tri_array,t_array,[points])
+
+    stage     = Stage.get_values(location='centroids',indices=[n_point])
+    xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])
+    ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])
+    #print '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum)
+    file.write( '%10.6f   %10.6f  %10.6f   %10.6f\n'%(t,stage,xmomentum,ymomentum) )
+
+    #time_array.append(t)
+    #stage_array.append(stage)
+
+    #if interactive_visualisation:
+    #    vis.update()
+"""
+file.close()
+print 'That took %.2f seconds' %(time.time()-t0)
+from pylab import *
+#ion()
+hold(False)
+plot(time_array, stage_array, 'r.-')
+#title('Gauge %s' %name)
+xlabel('time (s)')
+ylabel('stage (m)')
+#legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
+#savefig(name, dpi = 300)
+
+#raw_input('Next')
+show()
+"""
+
+#if interactive_visualisation:
+#    vis.evolveFinished()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/plot_results_cross_section.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/plot_results_cross_section.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/plot_results_cross_section.py	(revision 8734)
@@ -0,0 +1,58 @@
+"""
+    Quick plot of the dam break outputs
+    Note: v = p2_st.y[?], where ? is the cell position of the origin.
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from math import sqrt, pi,sin,cos
+from numpy import ones
+from analytical_paraboloid_basin import analytic_sol
+
+# Get the sww file
+p_st = util.get_output('paraboloid.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[11248]#[4898]
+v2=(p2_st.y==v)
+time_level = 50
+w2,u2,h2 = analytic_sol(p2_st.x[v2], v2, p2_st.time[time_level])
+
+
+#Plot stages
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[time_level,v2],'b.', label='numerical stage')
+pyplot.plot(p2_st.x[v2], p2_st.elev[v2],'k-', label='bed elevation')
+pyplot.plot(p2_st.x[v2], w2,'r-', label='analytical stage')
+pyplot.title('Stage at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.savefig('cross_section_stage.png')
+
+
+#Plot xmomentum
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[time_level,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], u2*h2,'r-', label='analytical')
+pyplot.title('Xmomentum at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('cross_section_xmom.png')
+
+
+#Plot velocities
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[time_level,v2], 'b.', label='numerical')
+pyplot.plot(p2_st.x[v2], u2,'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('cross_section_xvel.png')
+
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/plot_results_origin_wrt_time.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/plot_results_origin_wrt_time.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/plot_results_origin_wrt_time.py	(revision 8734)
@@ -0,0 +1,62 @@
+"""View results of runup.py
+"""
+#---------------
+# Import Modules
+#---------------
+#import anuga
+#import struct
+import numpy
+from matplotlib import pyplot as pyplot
+from anuga.utilities import plot_utils as util
+from numpy import sqrt, sin, cos, pi
+
+# Parameters
+D0 = 1000.
+L = 2500.
+R0 = 2000.
+g = 9.81
+
+A = (L**4 - R0**4)/(L**4 + R0**4)
+omega = 2./L*sqrt(2.*g*D0)
+T = pi/omega
+index_origin = 11248#4898
+
+
+# Get the sww file
+p = util.get_output('paraboloid.sww')
+p2=util.get_centroids(p)
+
+
+# Free surface over time at the origin of spatial domain.
+w = D0*((sqrt(1.-A*A))/(1.-A*cos(omega*p2.time)) -1.)
+
+#Plot stages
+pyplot.clf()
+pyplot.plot(p2.time,w, color='blue', label='analytical')
+pyplot.plot(p2.time,p2.stage[:,index_origin], color='green', label='numerical')
+pyplot.legend(loc='best')
+pyplot.title('Stage at the centre of the paraboloid basin over time')
+pyplot.xlabel('Time')
+pyplot.ylabel('Stage')
+pyplot.savefig('Stage_origin.png')
+
+#Plot xmomentums
+pyplot.clf()
+pyplot.plot(p2.time,0.0*p2.time, color='blue', label='analytical') # At the origin of space, w=h as always wet.
+pyplot.plot(p2.time,p2.xmom[:,index_origin], color='green', label='numerical')
+pyplot.legend(loc='best')
+pyplot.title('Xmomentum at the centre of the paraboloid basin over time')
+pyplot.xlabel('Time')
+pyplot.ylabel('Xmomentum')
+pyplot.savefig('Xmom_origin.png')
+
+#Plot velocities
+pyplot.clf()
+pyplot.plot(p2.time,0.0*p2.time, color='blue', label='analytical')
+pyplot.plot(p2.time,p2.xvel[:,index_origin], color='green', label='numerical')
+pyplot.legend(loc='best')
+pyplot.title('Xvelocity at the centre of the paraboloid basin over time')
+pyplot.xlabel('Time')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('Xvel_origin.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('run_paraboloid.py')
+    run_validation_script('plot_results_cross_section.py')
+    run_validation_script('plot_results_origin_wrt_time.py')    
+
+def clean():
+    autoclean()
+
+main()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/paraboloid_basin/results.tex	(revision 8734)
@@ -0,0 +1,60 @@
+\section{Oscillation on a paraboloid basin}
+This simulates flow oscillations in a paraboloid basin. The analytical solution is described by Thacker (REF), and is periodic. At any instant in time, the free surface elevation is paraboloid, and the velocity is linear. The scenario includes regular wetting and drying, as the flow oscillates back and forth in the basin. As well as testing the ability of the code to do wetting and drying, it will highlight any numerical energy loss or gain, manifest as an increase or decrease in the magnitude of the flow oscillations over long time periods, compared with the analytical solution. 
+
+\subsection{Results}
+
+Figures~\ref{fig:cs_stage}--\ref{fig:cs_xvel} show a good agreement between analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{cross_section_stage.png}
+\caption{Stage on a cross section of the basin at time $t=$ ?}
+\label{fig:cs_stage}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{cross_section_xmom.png}
+\caption{Xmomentum on a cross section of the basin at time $t=$ ?}
+\label{fig:cs_xmom}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{cross_section_xvel.png}
+\caption{Xvelocity on a cross section of the basin at time $t=$ ?}
+\label{fig:cs_xvel}
+\end{center}
+\end{figure}
+
+
+There should be a good agreement with the analytical solution, although as time goes on, some deviations may appear.  See Figures~\ref{fig:w_centre}--\ref{fig:u_centre}.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{Stage_centre.png}
+\caption{Stage over time in the centre of the paraboloid basin.}
+\label{fig:w_centre}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{Xmom_centre.png}
+\caption{Xmomentum over time in the centre of the paraboloid basin.}
+\label{fig:p_centre}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{Xvel_centre.png}
+\caption{Xvelocity over time in the centre of the paraboloid basin.}
+\label{fig:u_centre}
+\end{center}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/plot_results.py	(revision 8734)
@@ -0,0 +1,95 @@
+"""
+View results of run_channel.py
+
+Simple water flow example using ANUGA: Water flowing down a channel.
+It was called "steep_slope" in an old validation test.
+"""
+#---------------
+# Import Modules
+#---------------
+import time
+import anuga
+import numpy
+import scipy
+from matplotlib import pyplot as pyplot
+from anuga.utilities import plot_utils as util
+#--------------
+# Get variables
+#--------------
+p=util.get_output('channel.sww', 0.001)
+p2=util.get_centroids(p,velocity_extrapolation=True)
+
+#------------------
+# Select line
+#------------------
+#v=(y==2.5)
+v=(p2.y==p2.y[3])
+
+#-------------------------------
+# Define variables of case study
+#-------------------------------
+mann=0.03 # Manning's coef
+bedslope=-0.1
+fluxin=20./100. #The momentum flux at the upstream boundary ( = discharge / width)
+
+#---------------------------------------------
+# Analytical solution for steady, uniform flow
+#---------------------------------------------
+# This comes from solving 
+# 1) Friction_slope=bedslope , and
+# 2) flux = depth*velocity = flux_in
+uana= ( mann**(-2.)*abs(bedslope)*fluxin**(4./3.) )**(3./10.) # Velocity
+dana= fluxin/uana # Depth
+
+#--------------------
+# Make plot animation
+#--------------------
+pyplot.clf()
+pyplot.ion()
+
+line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,( (p2.stage[:,v]-p2.elev[:,v]).max(),(p2.stage[:,v]-p2.elev[v]).min() ) )
+line.set_label('numerical')
+pyplot.plot( (0,100),(dana,dana), 'r',label='analytical' )
+pyplot.legend()
+#pyplot.plot(x[v],elev[v],'green')
+for i in range(p2.xmom.shape[0]):
+    line.set_xdata(p2.x[v])
+    line.set_ydata(p2.stage[i,v]-p2.elev[v])
+    pyplot.draw()
+    pyplot.title('Flow depth: should be converging to steady uniform flow ' )
+    pyplot.xlabel('Xposition')
+    pyplot.ylabel('Depth')
+pyplot.ylim([0.085,0.095])
+pyplot.legend(loc='best')
+pyplot.title('Flow depth at 400s')# -- there should be a central region with steady uniform flow ' )
+pyplot.savefig('final_depth.png')
+
+#--------------------------------------------
+# Compare velocity with analytical solution
+#--------------------------------------------
+# Find an x value close to x==50
+tmp=(abs(p2.x-50.)).argmin()
+v=(p2.x==p2.x[tmp])
+
+pyplot.clf()
+pyplot.plot(p2.y[v],p2.xvel[100,v], 'o', label='numerical')
+pyplot.plot((0,100),(uana,uana),label='analytical')
+pyplot.ylim([1.0,3.0])
+pyplot.xlabel('Yposition along the line x=50')
+pyplot.ylabel('Xvelocity m/s')
+pyplot.title('Final Xvelocity around the line x=50.')
+pyplot.legend(loc='best')
+pyplot.savefig('x_velocity.png')
+
+pyplot.clf()
+pyplot.plot(p2.y[v],p2.yvel[100,v],'o', label='numerical')
+pyplot.plot((0,100),(0.0, 0.0),label='analytical')
+pyplot.xlabel('Yposition along the line x=50')
+pyplot.ylabel('Yvelocity')
+pyplot.title('Final Yvelocity around the line x=50.')
+pyplot.legend(loc='best')
+pyplot.savefig('y_velocity.png')
+
+
+
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/produce_results.py	(revision 8734)
@@ -0,0 +1,21 @@
+"""
+Simple water flow example using ANUGA: Water flowing down a channel.
+It was called "steep_slope" in an old validation test.
+"""
+#--------------------------------
+# import modules
+#--------------------------------
+
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('run_channel.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/results.tex	(revision 8734)
@@ -0,0 +1,39 @@
+%Note for this test case:
+%Simple water flow example using ANUGA: Water flowing down a channel.
+%It was called "steep_slope" in an old validation test.
+
+
+\section{Shallow flow down a mild slope}
+This case simulates very shallow flow running down a mild slope topography. It represents an idealisation of the rainfall-runoff problem, which will often involve very shallow flows down such a topography. This case has an analytical solution (the steady-uniform solution, bed slope = friction slope).   
+
+\subsection{Results}
+Figures~\ref{fig:depthdownchan} shows the steady state depth in the downstream direction. There should be a good agreement with the analytical solution, at least away from the boundaries.  
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.8\textwidth]{final_depth.png}
+\caption{Depth in the downstream direction}
+\label{fig:depthdownchan}
+\end{center}
+\end{figure}
+
+Figures~\ref{fig:xvelscrosschan} and~\ref{fig:yvelscroschan} show the steady state $x$- and $y$-velocities, along a slice in the cross slope direction (near $x=50$). The $x$-velocities should agree well with the analytical solution, and the $y$-velocities should be zero.  
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.8\textwidth]{x_velocity.png}
+\caption{$x$-velocity along the cross-section $x=50$ (i.e. a cross-section with constant bed elevation)}
+\label{fig:xvelscrosschan}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.8\textwidth]{y_velocity.png}
+\caption{$y$-velocity along the cross-section $x=50$ (i.e. a cross-section with constant bed elevation)}
+\label{fig:yvelscroschan}
+\end{center}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/run_channel.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/run_channel.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/rundown_mild_slope/run_channel.py	(revision 8734)
@@ -0,0 +1,66 @@
+"""
+Simple water flow example using ANUGA: Water flowing down a channel.
+It was called "steep_slope" in an old validation test.
+"""
+import sys
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import anuga
+from anuga import rectangular_cross as rectangular_cross
+from anuga.structures.inlet_operator import Inlet_operator
+from anuga import Domain
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+points, vertices, boundary = rectangular_cross(50, 50, len1=100.0, len2=100.0)
+domain = Domain(points, vertices, boundary) # Create domain
+domain.set_name('channel') # Output name
+
+#------------------------------------------------------------------------------
+# Setup Algorithm, either using command line arguments
+# or override manually yourself
+#------------------------------------------------------------------------------
+from anuga.utilities.argparsing import parse_standard_args
+alg, cfl = parse_standard_args()
+alg = '2_0'
+domain.set_flow_algorithm(alg)
+domain.set_CFL(cfl)
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+Qin=20.
+fluxin=Qin/100. #The momentum flux at the upstream boundary ( = discharge / width)
+uana= 2.15843634571 # analytical Xvelocity
+dana= 0.0926596702273 # analytical water depth
+
+def topography(x, y):
+	return -x/10. # linear bed slope
+
+def init_stage(x,y):
+    stg= -x/10.+0.01 # Constant depth: 1 cm.
+    return stg
+#line0=[ [0.,0.], [0., 100.] ]
+#Uin=[uana, 0.0]
+#Inlet_operator(domain, line0, Q=Qin, velocity=Uin)
+
+domain.set_quantity('elevation', topography) # Use function for elevation
+domain.set_quantity('friction', 0.03) # Constant friction
+domain.set_quantity('stage', init_stage)
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions
+#------------------------------------------------------------------------------
+Bt = anuga.Transmissive_boundary(domain)
+BdIN = anuga.Dirichlet_boundary([dana, fluxin, 0.0])
+#BdOUT = anuga.Dirichlet_boundary([dana-10., fluxin, 0.0])
+Br = anuga.Reflective_boundary(domain) # Solid reflective wall
+domain.set_boundary({'left': BdIN, 'right': Bt, 'top': Br, 'bottom': Br})
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=2.0, finaltime=200.0):
+	print domain.timestepping_statistics()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/plot_runup.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/plot_runup.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/plot_runup.py	(revision 8734)
@@ -0,0 +1,87 @@
+"""View results of runup.py
+"""
+#---------------
+# Import Modules
+#---------------
+import scipy
+import anuga
+from anuga.utilities import plot_utils as util
+from matplotlib import pyplot as pyplot
+from numpy import zeros
+
+p2=util.get_output('runup.sww', minimum_allowed_height=1.0e-03)
+p=util.get_centroids(p2, velocity_extrapolation=True)
+
+#------------------
+# Select line
+#------------------
+py_central=p.y[scipy.argmin(abs(p.y-0.5))]
+v=(p.y==p.y[py_central])
+
+#--------------------
+# Make plot animation
+#--------------------
+##pyplot.close() #If the plot is open, there will be problems
+##pyplot.ion()
+##if False:
+##    line, = pyplot.plot( (p.x[v].min(),p.x[v].max()) ,(p.xvel[:,v].min(),p.xvel[:,v].max() ) )
+##    for i in range(p.xmom.shape[0]):
+##        line.set_xdata(p.x[v])
+##        line.set_ydata(p.xvel[i,v])
+##        pyplot.draw()
+##        pyplot.plot( (0,1),(0,0), 'r' )
+##        pyplot.title(str(i)+'/200') # : velocity does not converge to zero' )
+##        pyplot.xlabel('x')
+##        pyplot.ylabel('Velocity (m/s)')
+##
+##    pyplot.savefig('runup_x_velocities.png')
+
+def analytic_sol(x):
+    y = zeros(len(x))
+    for i in range(len(x)):
+        if x[i] < 0.2:
+            y[i] = -0.5*x[i]
+        else:
+            y[i] = -0.1
+    return y
+
+#------------------------------------------------
+# Maximum y velocities -- occurs in output step 3
+#------------------------------------------------
+pyplot.clf()
+pyplot.plot(p.x[v],p.stage[5,v],'o',label='numerical stage')
+pyplot.plot(p.x[v],p.elev[v],'k-',label='bed elevation')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.title('Free surface and bed at y = 0.5, time = 1.0 second')
+pyplot.legend(loc='best')
+pyplot.savefig('stage_1s.png')
+
+pyplot.clf()
+pyplot.plot(p.x[v],p.xvel[5,v],'-o', label='numerical')
+pyplot.title('Xvelocity at y = 0.5, time = 1.0 second')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.legend(loc='best')
+pyplot.savefig('xvel_1s.png')
+
+
+
+pyplot.clf()
+pyplot.plot(p.x[v],p.stage[150,v],'o-',label='numerical stage')
+pyplot.plot(p.x[v], analytic_sol(p.x[v]), label='analytical stage')
+pyplot.plot(p.x[v],p.elev[v],'k-',label='bed elevation')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Stage')
+pyplot.title('Free surface and bed at y = 0.5, time = 30.0 second')
+pyplot.legend(loc='best')
+pyplot.savefig('stage_30s.png')
+
+pyplot.clf()
+pyplot.plot(p.x[v],p.xvel[150,v],'-o', label='numerical')
+pyplot.plot(p.x[v], zeros(len(p.x[v])), label='analytical')
+pyplot.legend(loc='best')
+pyplot.title('Xvelocity at y = 0.5, time = 30.0 second')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_30s.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/produce_results.py	(revision 8734)
@@ -0,0 +1,17 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('runup.py')
+    run_validation_script('plot_runup.py')
+
+def clean():
+    autoclean()
+
+main()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/results.tex	(revision 8734)
@@ -0,0 +1,41 @@
+\section{Simple wave runup}
+This scenario simulates a wave flowing up a planar beach. Following the initial wave runup, eventually the water elevation should become constant, and the velocities should approach zero. 
+
+\subsection{Results}
+Figure~\ref{fig:stage_1s} shows the water surface at time 1s (in the cross-shore direction). It is not constant as the water runup at this time. Figure~\ref{fig:xvel_1s} shows the corresponding $x$-velocity during the wave runup. The velocities should be free from major spikes.
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_1s.png}
+\caption{Water surface during the wave runup at time $t=1.0$\,.}
+\label{fig:stage_1s}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_1s.png}
+\caption{Xvelocity during the wave runup at time $t=1.0$\,.}
+\label{fig:xvel_1s}
+\end{center}
+\end{figure}
+
+
+
+Figure~\ref{fig:stage_30s} shows the water surface at time 30s (in the cross-shore direction). It should be nearly constant (= -0.1m) in the wet portions of the domain. Figure~\ref{fig:xvel_30s} shows the corresponding velocity at time 30s. It should be nearly zero (e.g. $<<$ 1 mm/s). This case has been used to illustrate wet-dry artefacts in some versions of \anuga.
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_30s.png}
+\caption{Water surface at time 30s after the wave runup. It should be nearly constant in wet parts of the domain.}
+\label{fig:stage_30s}
+\end{center}
+\end{figure}
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_30s.png}
+\caption{Xvelocity at time 30s after the wave runup. It should be nearly zero.}
+\label{fig:xvel_30s}
+\end{center}
+\end{figure}
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/runup.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/runup.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_beach/runup.py	(revision 8734)
@@ -0,0 +1,68 @@
+"""
+Runup example from the manual, slightly modified
+"""
+#---------
+#Import Modules
+#--------
+import anuga
+import numpy
+from math import sin, pi, exp
+from anuga import Domain
+
+#---------
+#Setup computational domain
+#---------
+points, vertices, boundary = anuga.rectangular_cross(100,3)
+domain=Domain(points,vertices,boundary)    # Create Domain
+domain.set_name('runup')                   # Output to file runup.sww
+domain.set_datadir('.')                    # Use current folder
+domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
+
+#------------------------------------------------------------------------------
+# 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)
+
+#------------------
+# Define topography
+#------------------
+def topography(x,y):
+        return -x/2 #Linear bed slope
+
+def stagefun(x,y):
+    return -0.2 #Stage
+
+domain.set_quantity('elevation',topography)     # Use function for elevation
+domain.get_quantity('elevation').smooth_vertex_values() # Steve's fix -- without this, substantial artificial velcities are generated everywhere in the domain. With this fix, there are artificial velocities near the coast, but not elsewhere.
+domain.set_quantity('friction',0.0)             # Constant friction
+domain.set_quantity('stage', stagefun)          # Constant negative initial stage
+
+#--------------------------
+# Setup boundary conditions
+#--------------------------
+Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
+Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
+Bd=anuga.Dirichlet_boundary([-0.2, 0., 0.])     # Constant boundary values -- not used in this example
+Bw=anuga.Time_boundary(domain=domain,
+	function=lambda t: [(0.0*sin(t*2*pi)-0.1)*exp(-t)-0.1,  0.0,  0.0]) # Time varying boundary -- get rid of the 0.0 to do a runup.
+
+#----------------------------------------------
+# Associate boundary tags with boundary objects
+#----------------------------------------------
+domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom':Br})
+
+#------------------------------
+#Evolve the system through time
+#------------------------------
+#xwrite=open("xvel.out","wb")
+#ywrite=open("yvel.out","wb")
+## Set print options to be compatible with file writing via the 'print' statement 
+#numpy.set_printoptions(threshold=numpy.nan, linewidth=numpy.nan)
+
+for t in domain.evolve(yieldstep=0.2,finaltime=30.0):
+    print domain.timestepping_statistics()
+print 'Finished'
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/plot_results.py	(revision 8734)
@@ -0,0 +1,78 @@
+"""View results of runup_sinusoid.py
+"""
+#---------------
+# Import Modules
+#---------------
+from anuga.utilities import plot_utils as util
+from matplotlib import pyplot as pyplot
+import numpy
+
+p1=util.get_output('runup_sinusoid.sww', 0.001)
+#p1=util.get_output('runup_sinusoid.sww')
+p2=util.get_centroids(p1, velocity_extrapolation=True)
+#p2=util.get_centroids(p1)
+
+#------------------
+# Select line
+#------------------
+#v=(p.y==0.5)
+v=(p2.y==p2.y[80])
+
+#--------------------
+# Make plot animation
+#--------------------
+pyplot.close() #If the plot is open, there will be problems
+##pyplot.ion()
+##if True:
+##    line, = pyplot.plot( (p2.x[v].min(),p2.x[v].max()) ,(p2.xvel[:,v].min(),p2.xvel[:,v].max() ) )
+##    for i in range(p2.xmom.shape[0]):
+##        line.set_xdata(p2.x[v])
+##        line.set_ydata(p2.xvel[i,v])
+##        pyplot.draw()
+##        pyplot.plot( (0,1),(0,0), 'r' )
+##        pyplot.title(str(i)+'/40') # : velocity does not converge to zero' )
+##        pyplot.xlabel('Xposition')
+##        pyplot.ylabel('Xvelocity')
+##    pyplot.savefig('runup_x_velocities.png')
+
+# Plot vertex values
+pyplot.clf()
+t1=numpy.argmin(abs(p1.time-1.0))
+t2=len(p1.time)-1
+pyplot.scatter(p1.x,p1.y,c=p1.elev,edgecolors='none', s=25)
+pyplot.colorbar()
+pyplot.quiver(p1.x,p1.y,p1.xvel[t1,:],p1.yvel[t1,:])
+pyplot.title('The maximum VERTEX speed is '+ str(p1.vel[t1,:].max()) + ' m/s at time '+ str(p1.time[t1])+' s')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Yposition')
+pyplot.savefig('vel_t1_vertex.png')
+
+# Plot vertex values
+pyplot.clf()
+pyplot.scatter(p2.x,p2.y,c=p2.elev,edgecolors='none', s=25)
+pyplot.colorbar()
+pyplot.quiver(p2.x,p2.y,p2.xvel[t1,:],p2.yvel[t1,:])
+pyplot.title('The maximum CENTROID speed is '+ str(p2.vel[t1,:].max()) + ' m/s at time ' + str(p1.time[t1]) + ' s')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Yposition')
+pyplot.savefig('vel_t1_centroid.png')
+
+# Plot vertex values
+pyplot.clf()
+pyplot.scatter(p1.x,p1.y,c=p1.elev,edgecolors='none', s=25)
+pyplot.colorbar()
+pyplot.quiver(p1.x,p1.y,p1.xvel[t2,:],p1.yvel[t2,:])
+pyplot.title('The maximum VERTEX speed is '+ str(p1.vel[t2,:].max()) + ' m/s at time ' + str(p1.time[t2]) +' s')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Yposition')
+pyplot.savefig('vel_t2_vertex.png')
+
+# Plot vertex values
+pyplot.clf()
+pyplot.scatter(p2.x,p2.y,c=p2.elev,edgecolors='none', s=25)
+pyplot.colorbar()
+pyplot.quiver(p2.x,p2.y,p2.xvel[t2,:],p2.yvel[t2,:])
+pyplot.title('The maximum CENTROID speed is '+ str(p2.vel[t2,:].max()) + ' m/s at time ' + str(p1.time[t2])+ ' s')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Yposition')
+pyplot.savefig('vel_t2_centroid.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('runup_sinusoid.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/results.tex	(revision 8734)
@@ -0,0 +1,25 @@
+\section{Wave runup over sinusoidal ridges}
+This scenario simulates a wave flowing up a beach with transverse sinusoidal ridges. Following the initial wave runup, eventually the water elevation should become constant, and the velocities should approach zero. 
+
+\subsection{Results}
+Figure~\ref{fig:vel_t1_centroid} shows the centroid velocities during the wave runup. The flow should be concentrating in the channels near the shore, and be free from major spikes.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{vel_t1_centroid.png}
+\caption{Velocity during the wave runup. Point color corresponds to the bed elevation}
+\label{fig:vel_t1_centroid}
+\end{center}
+\end{figure}
+
+Figure~\ref{fig:vel_t2_centroid} shows the velocities profile at time 40~s. They should be nearly zero (e.g. O($10^{-3}$) m$/s$). This case has been used to illustrate wet-dry artefacts in some versions of \anuga.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{vel_t2_centroid.png}
+\caption{Velocity at time 20s after the wave runup. The flow speed should be nearly zero. Point color corresponds to the bed elevation.}
+\label{fig:vel_t2_centroid}
+\end{center}
+\end{figure}
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/runup_sinusoid.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/runup_sinusoid.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/runup_on_sinusoid_beach/runup_sinusoid.py	(revision 8734)
@@ -0,0 +1,74 @@
+"""
+Runup example from the manual, slightly modified
+"""
+#---------
+#Import Modules
+#--------
+import anuga
+import numpy
+from math import sin, pi, exp
+from anuga import Domain
+
+#---------
+#Setup computational domain
+#---------
+domain = anuga.rectangular_cross_domain(40,40, len1=1., len2=1.)
+
+
+#------------------------------------------------------------------------------
+# 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)
+domain.set_name('runup_sinusoid')                         # Output to file runup.sww
+domain.set_datadir('.')                          # Use current folder
+domain.set_quantities_to_be_stored({'stage': 2, 'xmomentum': 2, 'ymomentum': 2, 'elevation': 1})
+domain.set_minimum_allowed_height(0.01)
+
+
+#------------------
+# Define topography
+#------------------
+scale_me=1.0
+def topography(x,y):
+	return (-x/2.0 +0.05*numpy.sin((x+y)*50.0))*scale_me
+
+def stagefun(x,y):
+    stge=-0.2*scale_me
+    return stge
+
+domain.set_quantity('elevation',topography)     # Use function for elevation
+domain.get_quantity('elevation').smooth_vertex_values() 
+domain.set_quantity('friction',0.0)             # Constant friction
+domain.set_quantity('stage', stagefun)          # Constant negative initial stage
+domain.get_quantity('stage').smooth_vertex_values()
+
+
+#--------------------------
+# Setup boundary conditions
+#--------------------------
+Br=anuga.Reflective_boundary(domain)            # Solid reflective wall
+Bt=anuga.Transmissive_boundary(domain)          # Continue all values of boundary -- not used in this example
+Bd=anuga.Dirichlet_boundary([-0.1*scale_me,0.,0.])       # Constant boundary values -- not used in this example
+
+#----------------------------------------------
+# Associate boundary tags with boundary objects
+#----------------------------------------------
+domain.set_boundary({'left': Br, 'right': Bd, 'top': Br, 'bottom':Br})
+
+#------------------------------
+#Evolve the system through time
+#------------------------------
+for t in domain.evolve(yieldstep=1.0,finaltime=40.0):
+    print domain.timestepping_statistics()
+    xx = domain.quantities['xmomentum'].centroid_values
+    yy = domain.quantities['ymomentum'].centroid_values
+    dd = domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values 
+    dd = (dd)*(dd>1.0e-03)+1.0e-06
+    vv = ( (xx/dd)**2 + (yy/dd)**2 )**0.5
+    vv = vv*(dd>1.0e-03)
+    print 'Peak velocity is: ', vv.max(), vv.argmax()
+print 'Finished'
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/analytical_subcritical.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/analytical_subcritical.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/analytical_subcritical.py	(revision 8734)
@@ -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
+
+qA = 4.42  # This is the imposed momentum
+hx = 2.0   # This is the water height downstream
+g  = 9.81  # Accelleration due to gravity
+
+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/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/numerical_subcritical.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/numerical_subcritical.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/numerical_subcritical.py	(revision 8734)
@@ -0,0 +1,110 @@
+"""
+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 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'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 0.1
+dy = dx
+L = 25.
+W = 3*dx
+
+# 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) 
+
+#------------------------------------------------------------------------------
+# 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
+#------------------------------------------------------------------------------
+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)
+
+#-----------------------------------------------------------------------------
+# 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})
+
+
+#===============================================================================
+##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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 1.0, finaltime = 300.):
+    #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/subcritical_over_bump/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/plot_results.py	(revision 8734)
@@ -0,0 +1,56 @@
+"""
+    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
+
+p_st = util.get_output('subcritical.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+h,z = analytic_sol(p2_st.x[v2])
+
+#Plot the stages##############################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[300,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 an instant in time')
+##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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[300,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 an instant in time')
+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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[300,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 4.42/h,'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_subcritical.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/subcritical_over_bump/results.tex	(revision 8734)
@@ -0,0 +1,35 @@
+
+\section{Flow over a bump: subcritical flow}
+
+This is a subcritical flow over a bump. No shock occurs.
+
+\subsection{Results}
+
+
+We should see excellent agreement between the analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/analytical_with_shock.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/analytical_with_shock.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/analytical_with_shock.py	(revision 8734)
@@ -0,0 +1,79 @@
+"""
+Transcritical flow over a bump with a shock.
+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, array, sqrt
+from scipy.optimize import fsolve
+from pylab import plot,show,ylim
+
+
+qA  = 0.18  # This is the imposed momentum
+hx  = 0.33  # This is the water height downstream
+g   = 9.81  # Accelleration due to gravity
+
+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)
+    zM= max(z)
+
+    def find_hM(hM): #to find the water height at the maxima of the bump
+        return hx**3 + (-qA**2/(2*g*hM**2)-hM-zM)*hx**2 + qA**2/(2*g)
+    hM = fsolve(find_hM, 0.5)
+
+    def find_h(h): #to find the water height at every spatial point after hM is found
+        return h**3 + (zb-qA**2/(2*g*hM**2)-hM-zM)*h**2 + qA**2/(2*g)
+
+    def find_hL(h): #to find the water height at every spatial point after hM is found
+        return h**3 + (zb-qA**2/(2*g*hx**2)-hx)*h**2 + qA**2/(2*g)
+
+    def find_shock(Q):
+        h1 = Q[0]
+        h2 = Q[1]
+        zs = Q[2]
+        F = zeros(3)
+        F[0] = h1**3 + (zs-qA**2/(2*g*hM**2)-hM-zM)*h1**2 + qA**2/(2*g)
+        F[1] = h2**3 + (zs-qA**2/(2*g*hx**2)-hx)*h2**2 + qA**2/(2*g)
+        F[2] = (qA**2)*(1.0/h1 - 1.0/h2) + 0.5*g*(h1**2 - h2**2)
+        return F
+
+    z_shock = fsolve(find_shock,array([0.2,0.3,0.005]))[2]
+
+    def find_zshock(x):
+        return 0.2 - 0.05*(x-10.0)**2 - z_shock
+    x_shock = fsolve(find_zshock,11.5)[0]
+
+    h=zeros(len(x))
+    for i in range(len(x)):
+        if x[i] <= 10.0:
+            zb   = z[i]
+            h[i] = fsolve(find_h,0.4)
+        elif x[i] <= x_shock:
+            zb   = z[i]
+            h[i] = fsolve(find_h,0.1)
+        elif x[i] <= 12.0:
+            zb   = z[i]
+            h[i] = fsolve(find_hL,0.3)
+        else:
+            h[i] = hx
+    return h,z
+        
+##N = 10001
+##L = 25.
+##x = linspace(0.0,L,N)
+##h,z = analytic_sol(x)
+##plot(x,h+z,'b.-', x,z,'k-')
+##ylim([0.0, 0.5])
+##show()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/numerical_transcritical.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/numerical_transcritical.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/numerical_transcritical.py	(revision 8734)
@@ -0,0 +1,111 @@
+"""
+Simple water flow example using ANUGA
+Transcritical flow over a bump with a shock
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+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 = 'transcritical_'+time
+output_dir = '.'
+output_file = 'transcritical'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 0.1
+dy = dx
+L = 25.
+W = 3*dx
+
+# 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) 
+
+#------------------------------------------------------------------------------
+# 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
+#------------------------------------------------------------------------------
+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 stage(x,y):
+    return 2.0*ones(len(x))
+domain.set_quantity('stage', stage)
+
+#-----------------------------------------------------------------------------
+# 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 
+BdL = anuga.Dirichlet_boundary([0.41373588752426715, 0.18, 0.]) # Constant boundary values
+BdR = anuga.Dirichlet_boundary([0.33, 0.18, 0.]) # Constant boundary values
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': BdL, 'right': BdR, '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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 1.0, finaltime = 300.):
+    #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/transcritical_with_shock/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/plot_results.py	(revision 8734)
@@ -0,0 +1,55 @@
+"""
+    Quick plot for outputs of the transcritical flow with a shock
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_with_shock import *
+from numpy import ones
+
+p_st = util.get_output('transcritical.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+h,z = analytic_sol(p2_st.x[v2])
+
+#Plot the stages##############################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[300,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 an instant in time')
+##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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[300,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 0.18*ones(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.savefig('xmom_plot.png')
+
+
+
+#Plot the velocities#########################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[300,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 0.18/h,'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_transcritical.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_with_shock/results.tex	(revision 8734)
@@ -0,0 +1,35 @@
+
+\section{Flow over a bump: transcritical flow with a shock}
+
+This is transcritical flow over a bump with a shock.
+
+\subsection{Results}
+
+
+We should see excellent agreement between the analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/analytical_without_shock.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/analytical_without_shock.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/analytical_without_shock.py	(revision 8734)
@@ -0,0 +1,54 @@
+"""
+Transcritical flow over a bump without a shock.
+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, show
+
+
+q0  = 1.53  # This is the imposed momentum
+h_d = 0.66  # This is the water height downstream
+g   = 9.81  # Accelleration due to gravity
+
+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)
+    zM= max(z)
+
+    def find_hM(hM): #to find the water height at the maxima of the bump
+        return h_d**3 + (-q0**2/(2*g*hM**2)-hM-zM)*h_d**2 + q0**2/(2*g)
+    hM = fsolve(find_hM, 0.5)
+
+    def find_h(h): #to find the water height at every spatial point after hM is found
+        return h**3 + (zb-q0**2/(2*g*hM**2)-hM-zM)*h**2 + q0**2/(2*g)
+    h = zeros(len(x))
+    for i in range(len(x)):
+        zb = z[i]
+        #h[i] = fsolve(find_h, 1.0)
+        if x[i] < 10:
+            h[i] = fsolve(find_h, 1.0)
+        else:
+            h[i] = fsolve(find_h, 0.4)
+    return h, z
+
+##N = 401
+##L = 25.
+##x = linspace(0.0,L,N)
+##h,z=analytic_sol(x)
+##plot(x,h+z, x,z)
+##plot(x, 1.53/h)
+##show()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/numerical_transcritical.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/numerical_transcritical.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/numerical_transcritical.py	(revision 8734)
@@ -0,0 +1,119 @@
+"""
+Simple water flow example using ANUGA:
+Transcritical flow without shock over a bump.
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+import sys
+import anuga
+from anuga import Domain as Domain
+from math import cos
+from numpy import zeros, ones, float
+from time import localtime, strftime, gmtime
+from anuga.operators.set_w_uh_vh_operators import Polygonal_set_w_uh_vh_operator
+
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#-------------------------------------------------------------------------------
+time = strftime('%Y%m%d_%H%M%S',localtime())
+
+#output_dir = 'transcritical_'+time
+output_dir = '.'
+output_file = 'transcritical'
+
+#anuga.copy_code_files(output_dir,__file__)
+#start_screen_catcher(output_dir+'_')
+
+
+#------------------------------------------------------------------------------
+# Setup domain
+#------------------------------------------------------------------------------
+dx = 0.1
+dy = dx
+L = 25.
+W = 3*dx
+
+BC_polygonL = [[0,0],[dx,0.0], [dx,W], [0,W]]
+BC_polygonR = [[L,0.],[L,W], [L-dx,W], [L-dx,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) 
+
+#------------------------------------------------------------------------------
+# 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
+#------------------------------------------------------------------------------
+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
+        else:
+            z_b[i] = 0.0
+    return z_b
+domain.set_quantity('elevation',elevation)
+domain.set_quantity('friction', 0.0)
+domain.set_quantity('xmomentum', 0.0)
+
+
+def stage(x,y):
+    return 0.66*ones(len(x))
+domain.set_quantity('stage', stage)
+
+#-----------------------------------------------------------------------------
+# 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 
+BdL = anuga.Dirichlet_boundary([1.0144468506259066, 1.53, 0.]) # Constant boundary values
+BdR = anuga.Dirichlet_boundary([0.4057809296474606, 1.53, 0.]) # Constant boundary values
+
+
+# Associate boundary tags with boundary objects
+domain.set_boundary({'left': BdL, 'right': BdR, 'top': Br, 'bottom': Br})
+
+w_uh_vhL= [1.0144468506259066, 1.53, 0.]
+w_uh_vhR= [0.4057809296474606, 1.53, 0.]
+Polygonal_set_w_uh_vh_operator(domain,w_uh_vhL,BC_polygonL)
+Polygonal_set_w_uh_vh_operator(domain,w_uh_vhR,BC_polygonR)
+
+#===============================================================================
+##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()
+#===============================================================================
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep = 1.0, finaltime = 300.):
+    #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/transcritical_without_shock/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/plot_results.py	(revision 8734)
@@ -0,0 +1,56 @@
+"""
+    Quick plot for the outputs of transcritical flow without shock
+
+"""
+import anuga.utilities.plot_utils as util
+from matplotlib import pyplot as pyplot
+from analytical_without_shock import *
+from numpy import ones
+
+p_st = util.get_output('transcritical.sww')
+p2_st=util.get_centroids(p_st)
+
+v = p2_st.y[10]
+v2=(p2_st.y==v)
+
+h,z = analytic_sol(p2_st.x[v2])
+
+#Plot the stages##############################################################
+pyplot.clf()
+pyplot.ion()
+pyplot.plot(p2_st.x[v2], p2_st.stage[300,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 an instant in time')
+##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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xmom[300,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 1.53*ones(len(p2_st.x[v2])),'r-', label='analytical')
+pyplot.title('Xmomentum at an instant in time')
+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.ion()
+pyplot.plot(p2_st.x[v2], p2_st.xvel[300,v2], 'b.-', label='numerical') # 0*T/6
+pyplot.plot(p2_st.x[v2], 1.53/h,'r-', label='analytical')
+pyplot.title('Xvelocity at an instant in time')
+pyplot.legend(loc='best')
+pyplot.xlabel('Xposition')
+pyplot.ylabel('Xvelocity')
+pyplot.savefig('xvel_plot.png')
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/produce_results.py	(revision 8734)
@@ -0,0 +1,18 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('numerical_transcritical.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
+
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/transcritical_without_shock/results.tex	(revision 8734)
@@ -0,0 +1,35 @@
+
+\section{Flow over a bump: transcritical flow without a shock}
+
+This is transcritical flow over a bump without a shock.
+
+\subsection{Results}
+
+
+We should see excellent agreement between the analytical and numerical solutions.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{stage_plot.png}
+\end{center}
+\caption{Stage results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xmom_plot.png}
+\end{center}
+\caption{Xmomentum results}
+\end{figure}
+
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{xvel_plot.png}
+\end{center}
+\caption{Xvelocity results}
+\end{figure}
+
+
+\endinput
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/channel_floodplain.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/channel_floodplain.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/channel_floodplain.py	(revision 8734)
@@ -0,0 +1,188 @@
+"""
+Simple water flow example using ANUGA
+Water flowing down a channel with a floodplain
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+# Import standard shallow water domain and standard boundaries.
+import anuga
+import numpy
+from anuga.structures.inlet_operator import Inlet_operator
+from anuga import *
+
+#------------------------------------------------------------------------------
+# Useful parameters for controlling this case
+#------------------------------------------------------------------------------
+floodplain_length = 1000.0 # Model domain length
+floodplain_width = 14.0    # Model domain width
+floodplain_slope = 1./300.
+chan_initial_depth = 0.65  # Initial depth of water in the channel
+chan_bankfull_depth = 1.0  # Bankfull depth of the channel
+chan_width = 10.0          # Bankfull width of the channel
+bankwidth = 2.             # Width of the bank regions -- note that these protrude into the channel
+man_n=0.03                 # Manning's n
+l0 = 1.000 # Length scale associated with triangle side length in channel (min_triangle area = 0.5*l0^2)
+
+assert chan_width < floodplain_width, \
+        'ERROR: Channel width is greater than floodplain width'
+
+assert bankwidth < chan_width/2., \
+        'ERROR: The bank width must be less than half the channel width'
+
+#------------------------------------------------------------------------------
+# Setup computational domain
+#------------------------------------------------------------------------------
+# Define boundary polygon -- in this case clockwise around the proposed boundary
+boundary_polygon = [ [0.,0.], 
+                     [0., floodplain_length], 
+                     [floodplain_width/2. - chan_width/2., floodplain_length], 
+                     [floodplain_width/2. + chan_width/2., floodplain_length], 
+                     [floodplain_width, floodplain_length], 
+                     [floodplain_width, 0.], 
+                     [floodplain_width/2. + chan_width/2., 0.], 
+                     [floodplain_width/2. - chan_width/2., 0.] 
+                     ]
+
+# Define channel polygon, where we can optionally refine the resolution. 
+# Note that we set this a distance l0 inside the boundary polygon, so the polygons
+# do not overlap. 
+channel_polygon = [ [floodplain_width/2. - chan_width/2., +l0],
+                    [floodplain_width/2. - chan_width/2., floodplain_length-l0],
+                    [floodplain_width/2. + chan_width/2., floodplain_length-l0],
+                    [floodplain_width/2. + chan_width/2., +l0]
+                    ]
+
+# Define domain with appropriate boundary conditions
+domain = create_domain_from_regions( boundary_polygon, 
+                                   boundary_tags={'left': [0],
+                                                  'top1': [1],
+                                                  'chan_out': [2],
+                                                  'top2': [3],
+                                                  'right': [4],
+                                                  'bottom1': [5],
+                                                  'chan_in': [6],
+                                                  'bottom2': [7] },
+                                   maximum_triangle_area = 0.5*l0*l0,
+                                   minimum_triangle_angle = 28.0,
+                                   mesh_filename = 'channel_floodplain.msh',
+                                   interior_regions = [ ],
+                                   #interior_regions = [\
+                                   #    [channel_polygon, 0.5*l0*l0] ],
+                                   use_cache=False,
+                                   verbose=True)
+
+domain.set_name('channel_floodplain') # Output name
+
+#------------------------------------------------------------------------------
+# 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
+#
+#------------------------------------------------------------------------------
+
+# Function for topography
+def elevation(x, y):
+    elev1= -y*floodplain_slope - chan_bankfull_depth*\
+            (x>(floodplain_width/2. - chan_width/2.))*\
+            (x<(floodplain_width/2. + chan_width/2.)) 
+    # Add banks
+    if(bankwidth>0.0):
+        leftbnk = floodplain_width/2. - chan_width/2.
+        rightbnk = floodplain_width/2. + chan_width/2.
+        # Left bank
+        elev2 = elev1 + (chan_bankfull_depth \
+                - chan_bankfull_depth/bankwidth*(x - leftbnk))*\
+                (x>leftbnk)*(x < leftbnk + bankwidth)
+        # Right bank
+        elev2 = elev2 + (chan_bankfull_depth \
+                + chan_bankfull_depth/bankwidth*(x - rightbnk))*\
+                (x>rightbnk-bankwidth)*(x < rightbnk)
+    if(bankwidth==0.0):
+        elev2 = elev1
+    return elev2
+
+
+#Function for stage
+def stage(x,y):
+    return -y*floodplain_slope -chan_bankfull_depth + chan_initial_depth 
+
+domain.set_quantity('elevation', elevation) # Use function for elevation
+domain.set_quantity('friction', man_n)      # Constant friction
+domain.set_quantity('stage', stage)         # Use function for stage
+
+# Define inlet operator 
+flow_in_yval=0.0
+if True:
+    line1 = [ [floodplain_width/2. - chan_width/2., flow_in_yval],\
+              [floodplain_width/2. + chan_width/2., flow_in_yval] ]
+    Qin = 0.5*(floodplain_slope*(chan_width*chan_initial_depth)**2.*man_n**(-2.)\
+            *chan_initial_depth**(4./3.) )**0.5
+    Inlet_operator(domain, line1, Qin)
+    print 'Discharge in = ', Qin 
+
+#------------------------------------------------------------------------------
+#
+# Setup boundary conditions
+#
+#------------------------------------------------------------------------------
+
+Br = anuga.Reflective_boundary(domain) # Solid reflective wall
+
+def outflow_stage_boundary(t):
+    return -floodplain_length*floodplain_slope \
+            + chan_initial_depth - chan_bankfull_depth
+
+# Note that the outflow boundary may be slightly incorrect for the trapezoidal channel case, 
+# or incorrect more generally if there are numerical problems. But, in the central regions of
+# the channel, this shouldn't prevent us reaching steady, uniform flow.
+Bout_tmss = anuga.shallow_water.boundaries.Transmissive_momentum_set_stage_boundary(domain, function = outflow_stage_boundary) 
+
+domain.set_boundary({'left': Br, 
+                     'right': Br, 
+                     'top1': Br, 
+                     'top2': Br, 
+                     'bottom1': Br, 
+                     'bottom2': Br, 
+                     'chan_out': Bout_tmss, 
+                     'chan_in': Br})
+
+# Set up file to record computations of discharge at several points.
+discharge_outfile=open('discharge_outputs.txt', 'w')
+discharge_outfile.write('Time (s)'+","+ 'Discharge@10' + ","+ 'Discharge@700'+","+ 'Discharge@1000' + "\n")
+
+#------------------------------------------------------------------------------
+#
+# Evolve system through time
+#
+#------------------------------------------------------------------------------
+for t in domain.evolve(yieldstep=2.0, finaltime=1200.0):
+    print domain.timestepping_statistics()
+    xx=domain.quantities['ymomentum'].centroid_values
+    dd=(domain.quantities['stage'].centroid_values - domain.quantities['elevation'].centroid_values)
+    dd=dd*(dd>0.)
+
+    tmp = xx/(dd+1.0e-06)*(dd>0.0)
+    print tmp.max(), tmp.argmax(), tmp.min(),  tmp.argmin()
+
+    # Compute flow through cross-section -- check that the inflow boundary condition is doing its job
+    # This also provides another useful steady-state check
+    if( numpy.floor(t/100.) == t/100. ):
+        print '#### COMPUTING FLOW THROUGH CROSS-SECTIONS########'
+        s0 = domain.get_flow_through_cross_section([[0., 10.0], [floodplain_width, 10.0]])
+        s1 = domain.get_flow_through_cross_section([[0., floodplain_length-300.0], [floodplain_width, floodplain_length-300.0]])
+        s2 = domain.get_flow_through_cross_section([[0., floodplain_length-1.0], [floodplain_width, floodplain_length-1.0]])
+        
+        print 'Cross sectional flows: ',s0, s1, s2 
+        print '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
+        discharge_outfile.write(str(t) + "," +str(s0) + ","+ str(s1) +"," + str(s2) + "\n")
+discharge_outfile.close()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/plot_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/plot_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/plot_results.py	(revision 8734)
@@ -0,0 +1,141 @@
+from anuga.utilities import plot_utils as util
+from matplotlib import pyplot as pyplot
+
+# Time-index to plot outputs from
+index=600
+p2 = util.get_output('channel_floodplain.sww')
+p=util.get_centroids(p2, velocity_extrapolation=True)
+v = (p.x>6.0)*(p.x<8.0)
+
+# Numerical results along a central channel 'slice'
+V1 = p.stage[index,v] - p.elev[v]
+V2 = p.yvel[index,v]
+V3 = p.xvel[index,v]
+
+##########################################################################
+# Analytical solution of steady uniform 2D flow in a trapezoidal channel.
+##########################################################################
+
+Qin=4.6932   # Inflow discharge
+slp=1./300.  # Floodplain slope (= water slope for steady uniform flow)
+man_n=0.03   # Manning's n
+Bcentral=6.0 #Flat bed width of the trapezoidal channel
+alpha=0.5    # Side slope of the trapezoidal banks
+
+k = (slp*(1./man_n)**2)**0.5 # At any point, the analytical solution says U = k*d^(2/3)
+
+# Function to calculate the discharge, given the channel centre depth dc, assuming
+# steady uniform flow
+def discharge_su(dc):
+    if(alpha>0.):
+        out = 2*k*( 3./(8.*alpha)*(dc)**(8./3.)) +Bcentral*k*(dc)**(5./3.)
+    else:
+        out = Bcentral*k*(dc)**(5./3.)    
+    return out 
+
+# Function that will be minimized to find the depth associated with discharge Qin
+def minme(dc):
+    q1 = discharge_su(dc)
+    return (q1-Qin)**2.
+
+# Minimise the function mimne, to find the centre depth.
+import scipy.optimize
+dc_analytical = scipy.optimize.fmin(minme, x0=1.0)[0]
+
+
+
+##################################
+# Plots
+##################################
+
+# Analytical solution has U*abs(U)*n^2 / D^(4./3.) = Sf = bed slope
+# Hence, the following two variables should be equal -- I have checked that x velocities are fairly small
+pyplot.clf()
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.y[v], V2*0+k*dc_analytical**(2./3.), 'o', label='analytical velocity')
+pyplot.plot(p.y[v], (V2**2)**0.5,'o', label='numerical velocity')
+pyplot.plot(p.y[v], V1, 'o',label='numerical depth')
+pyplot.plot(p.y[v], V1*0. + dc_analytical, 'o', label='analytical depth')
+pyplot.title('Mid channel numerical velocities and depths, vs analytical velocities and depths')
+pyplot.legend(loc='best')
+pyplot.xlabel('Down-channel distance (m)')
+pyplot.ylabel('Generic scale (m or m/s)')
+pyplot.savefig('fig1mid_channel.png')
+
+
+
+# Plot velocity over the cross-section
+pyplot.clf()
+v1 = (p.y<105.0)&(p.y>95.0)
+analytical_stage = min(p.elev[v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+
+temp0 = p.stage[index,v1]*0. + analytical_stage
+temp1 = (temp0) * (temp0 > p.elev[v1])
+temp2 = (p.elev[v1]) * (temp0 < p.elev[v1])
+Analytic_Stage = temp1 + temp2
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity')
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='numerical velocity')
+#pyplot.plot(p.x[v1],p.stage[index,v1]-p.elev[v1],'ko', label='numerical height')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='numerical stage')
+pyplot.plot(p.x[v1],Analytic_Stage,'o', label='analytical stage')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation')
+pyplot.legend(loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Downstream channel regions (95 to 105m)' +'\n')
+pyplot.xlabel('Cross-channel distance (m)')
+pyplot.ylabel('Generic scale (m or m/s)')
+pyplot.savefig('fig2upstream_channel.png')
+
+
+
+# Plot velocity over the cross-section
+pyplot.clf()
+v1 = (p.y<505.0)&(p.y>495.0)
+analytical_stage = min(p.elev[v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+
+temp0 = p.stage[index,v1]*0. + analytical_stage
+temp1 = (temp0) * (temp0 > p.elev[v1])
+temp2 = (p.elev[v1]) * (temp0 < p.elev[v1])
+Analytic_Stage = temp1 + temp2
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity')
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='numerical velocity')
+#pyplot.plot(p.x[v1],p.stage[index,v1]-p.elev[v1],'ko', label='numerical height')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='numerical stage')
+pyplot.plot(p.x[v1],Analytic_Stage,'o', label='analytical stage')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation')
+pyplot.legend(loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Central channel regions (495 to 505m)' +'\n')
+pyplot.xlabel('Cross-channel distance (m)')
+pyplot.ylabel('Generic scale (m or m/s)')
+pyplot.savefig('fig3central_channel.png') 
+
+
+
+# Plot velocity over the cross-section
+pyplot.clf()
+v1 = (p.y<905.0)&(p.y>895.0)
+analytical_stage = min(p.elev[v1]) + dc_analytical
+analytic_vel = ( (1./300.)*(analytical_stage-p.elev[v1])**(4./3.)*(1./0.03)**2.)**0.5
+analytic_vel = analytic_vel*(analytical_stage>p.elev[v1])
+
+temp0 = p.stage[index,v1]*0. + analytical_stage
+temp1 = (temp0) * (temp0 > p.elev[v1])
+temp2 = (p.elev[v1]) * (temp0 < p.elev[v1])
+Analytic_Stage = temp1 + temp2
+pyplot.figure(figsize=(12.,8.))
+pyplot.plot(p.x[v1], analytic_vel,'o', label='analytical velocity')
+pyplot.plot(p.x[v1], p.yvel[index,v1],'o', label='numerical velocity')
+#pyplot.plot(p.x[v1],p.stage[index,v1]-p.elev[v1],'ko', label='numerical height')
+pyplot.plot(p.x[v1],p.stage[index,v1],'o', label='numerical stage')
+pyplot.plot(p.x[v1],Analytic_Stage,'o', label='analytical stage')
+pyplot.plot(p.x[v1],p.elev[v1],'o', label='bed elevation')
+pyplot.legend(loc=10)
+pyplot.title('Velocity (analytical and numerical) and Stage:' + '\n' +'Downstream channel regions (895 to 905m)' +'\n')
+pyplot.xlabel('Cross-channel distance (m)')
+pyplot.ylabel('Generic scale (m or m/s)')
+pyplot.savefig('fig4downstream_channel.png')
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/produce_results.py
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/produce_results.py	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/produce_results.py	(revision 8734)
@@ -0,0 +1,16 @@
+#--------------------------------
+# import modules
+#--------------------------------
+from fabricate import *
+from validation_tests.utilities import run_validation_script
+
+# Setup the python scripts which produce the output for this
+# validation test
+def build():
+    run_validation_script('channel_floodplain.py')
+    run_validation_script('plot_results.py')
+
+def clean():
+    autoclean()
+
+main()
Index: /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/results.tex
===================================================================
--- /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/results.tex	(revision 8734)
+++ /trunk/anuga_core/source/anuga_validation_tests/Analytical_exact/trapezoidal_channel/results.tex	(revision 8734)
@@ -0,0 +1,38 @@
+\section{Trapezoidal channel steady uniform flow case}
+This scenario simulates steady, uniform flow in a trapezoidal channel. It tests the ability of the model to compute the correct flow depth, and the correct lateral distribution of flow velocity, and to conserve mass.
+
+\subsection{Results}
+Figure~\ref{fig:fig1mid_channel} shows the along-channel depth and velocity (in the deepest part of the channel). Both should be nearly constant away from the boundaries, and be close to the analytical solutions. 
+
+We do not expect perfect agreement, because the mesh is not very fine in this example (triangle side length of around 1m, just enough to resolve the banks). There will probably be some numerical diffusion in the cross-channel velocity profiles, which will in turn cause errors in the mid-channel velocity and free surface elevation. We deliberately choose to not use a finer mesh, because in realistic problems, it is often not possible to resolve all channels very well.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{fig1mid_channel.png}
+\caption{Depth and $y$-velocity in the middle of the channel, in the downstream direction.}
+\label{fig:fig1mid_channel}
+\end{center}
+\end{figure}
+
+Figure~\ref{fig:xsect_vels} shows the cross-channel velocity profiles at a number of cross-sections. Ideally it should agree with the analytical solution, however, this may be difficult due to numerical diffusion in the cross-channel direction. Irrespective, the velocity profile should be qualitatively correct -- highest velocities should be in the channel centre, with lower velocities towards the banks.
+
+\begin{figure}[h]
+\begin{center}
+\includegraphics[width=0.9\textwidth]{fig2upstream_channel.png}
+\includegraphics[width=0.9\textwidth]{fig3central_channel.png}
+\includegraphics[width=0.9\textwidth]{fig4downstream_channel.png}
+\caption{$y$-velocity distribution over a number of cross-sections.}
+\label{fig:xsect_vels}
+\end{center}
+\end{figure}
+
+Table~\ref{tab:trapztab} shows the discharge computed at a number of cross-sections in the channel, at a number of time-steps on the way to near steady-state. By the end of the simulation they should all be essentially the same. Large variations may suggest mass conservation errors (small variations are probably due to the interpolation that occurs in the 'compute\_flow\_through\_cross\_section' routine).
+
+\DTLloaddb{dischargeout}{Tests/Simple/trapezoidal_channel/discharge_outputs.txt}
+\begin{table}
+\caption{Discharge through cross-sections at a number of $x$-position, at different instants in time}
+\label{tab:trapztab}
+\DTLdisplaydb{dischargeout}
+\end{table}
+
+\endinput
