Index: /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/analytical_MacDonald.py
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/analytical_MacDonald.py	(revision 8677)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/analytical_MacDonald.py	(revision 8677)
@@ -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/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/numerical_MacDonald.py
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/numerical_MacDonald.py	(revision 8677)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/numerical_MacDonald.py	(revision 8677)
@@ -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/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/plot_results.py
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/plot_results.py	(revision 8677)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/plot_results.py	(revision 8677)
@@ -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/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/produce_results.py
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/produce_results.py	(revision 8677)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/produce_results.py	(revision 8677)
@@ -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/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/results.tex
===================================================================
--- /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/results.tex	(revision 8677)
+++ /trunk/anuga_core/validation_tests/Tests/Analytical_exact/MacDonald_short_channel/results.tex	(revision 8677)
@@ -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
