Index: /anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py
===================================================================
--- /anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py	(revision 3916)
+++ /anuga_validation/automated_validation_tests/okushiri_tank_validation/compare_timeseries_with_measures.py	(revision 3916)
@@ -0,0 +1,404 @@
+"""Verify that simulation produced by ANUGA compares to published
+validation timeseries ch5, ch7 and ch9 as well as the boundary timeseries.
+
+RMS norm is printed and plots are produced.
+"""
+
+from Numeric import allclose, argmin, argmax
+from Scientific.IO.NetCDF import NetCDFFile
+
+from anuga.abstract_2d_finite_volumes.util import file_function
+from anuga.utilities.numerical_tools import\
+     ensure_numeric, cov, get_machine_precision
+
+import project
+
+try:
+    from pylab import ion, hold, plot, title, legend, xlabel, ylabel, savefig
+except:
+    plotting = False
+else:
+    plotting = True
+
+#plotting = False
+
+#-------------------------
+# Basic data
+#-------------------------
+
+finaltime = 22.5
+timestep = 0.05
+
+gauge_locations = [[0.000, 1.696]] # Boundary gauge
+gauge_locations += [[4.521, 1.196],  [4.521, 1.696],  [4.521, 2.196]] #Ch 5-7-9
+gauge_names = ['Boundary', 'ch5', 'ch7', 'ch9']
+
+validation_data = {}
+for key in gauge_names:
+    validation_data[key] = []
+
+
+# Results from lwru2_variable_mesh.sww
+#
+#Validating Boundary
+#Covariance = 5.265850426352561254e-05
+#Accumulated difference = 9.188561500629542633e-04
+#Difference in maxima = 4.865498851439054029e-05
+#Difference in minima = 2.293007809784416290e-04
+#Timelag between maxima = 3.000000000000007105e-01
+#Timelag between minima = 3.000000000000007105e-01
+#Next
+#
+#Validating ch5
+#Covariance = 1.171913883449135149e-04
+#Accumulated difference = 3.439740513194261117e-03
+#Difference in maxima = 1.808007760316858448e-03
+#Difference in minima = 1.494336068527772621e-04
+#Timelag between maxima = 9.999999999999786837e-02
+#Timelag between minima = 8.499999999999996447e-01
+#Next#
+#
+#Validating ch7
+#Covariance = 1.130834789656423632e-04
+#Accumulated difference = 2.868914369295409543e-03
+#Difference in maxima = 5.464055727564115506e-04
+#Difference in minima = 4.515664278503157131e-03
+#Timelag between maxima = 0.000000000000000000e+00
+#Timelag between minima = 4.500000000000010658e-01
+#
+#Next
+#Validating ch9
+#Covariance = 1.257648236914514148e-04
+#Accumulated difference = 3.287722754612562841e-03
+#Difference in maxima = 1.856547605566957748e-03
+#Difference in minima = 3.406785970601119290e-03
+#Timelag between maxima = 0.000000000000000000e+00
+#Timelag between minima = 9.499999999999992895e-01
+
+
+
+
+
+# New validation mesh
+#Validating Boundary
+#Covariance = 5.294418463841652839e-05
+#Accumulated difference = 8.222031767747334978e-04
+#Difference in maxima = 2.074639250243537347e-05
+#Difference in minima = 2.594698240608837858e-04
+#Timelag between maxima = 2.500000000000000000e-01
+#Timelag between minima = 2.500000000000000000e-01
+#Next
+#
+#Validating ch5
+#Covariance = 1.169700827457795071e-04
+#Accumulated difference = 3.390769873961110153e-03
+#Difference in maxima = 1.443103996215826246e-03
+#Difference in minima = 6.360676810760659133e-04
+#Timelag between maxima = 4.999999999999715783e-02
+#Timelag between minima = 7.999999999999989342e-01
+#Next
+#
+#Validating ch7
+#Covariance = 1.131522443178145860e-04
+#Accumulated difference = 2.815084054773397005e-03
+#Difference in maxima = 8.380973795966278894e-05
+#Difference in minima = 4.604184345996111850e-03
+#Timelag between maxima = 0.000000000000000000e+00
+#Timelag between minima = 2.050000000000000711e+00
+#
+#Next
+#Validating ch9
+#Covariance = 1.251932626774203753e-04
+#Accumulated difference = 3.208757435748204462e-03
+#Difference in maxima = 2.081659579598611753e-03
+#Difference in minima = 3.428854197019200710e-03
+#Timelag between maxima = 0.000000000000000000e+00
+#Timelag between minima = 1.000000000000000000e+00
+
+
+expected_covariance = {'Boundary': 5.294418463841652839e-05, 
+                       'ch5': 1.169700827457795071e-04,
+                       'ch7': 1.131522443178145860e-04,
+                       'ch9': 1.251932626774203753e-04}
+
+expected_difference = {'Boundary': 8.222031767747334978e-04,
+                       'ch5': 3.390769873961110153e-03,
+                       'ch7': 2.815084054773397005e-03,
+                       'ch9': 3.208757435748204462e-03}
+
+expected_maximum_diff = {'Boundary': 2.074639250243537347e-05,
+                         'ch5': 1.443103996215826246e-03,
+                         'ch7': 8.380973795966278894e-05,
+                         'ch9': 2.081659579598611753e-03}
+
+expected_minimum_diff = {'Boundary': 2.594698240608837858e-04,
+                         'ch5': 6.360676810760659133e-04,
+                         'ch7': 4.604184345996111850e-03,
+                         'ch9': 3.428854197019200710e-03}
+
+expected_argmax_timelag = {'Boundary': 2.500000000000000000e-01,
+                           'ch5': 4.999999999999715783e-02,
+                           'ch7': 0.000000000000000000e+00,
+                           'ch9': 0.000000000000000000e+00}
+
+expected_argmin_timelag = {'Boundary': 2.500000000000000000e-01,
+                           'ch5': 7.999999999999989342e-01,
+                           'ch7': 2.050000000000000711e+00,
+                           'ch9': 1.000000000000000000e+00}
+
+
+
+
+# New validation mesh - smaller regions
+#Validating Boundary
+#Covariance = 5.269749049737429159e-05
+#Accumulated difference = 8.337114775060755880e-04
+#Difference in maxima = 7.166795837686146253e-05
+#Difference in minima = 1.335627610594206094e-04
+#Timelag between maxima = 3.000000000000007105e-01
+#Timelag between minima = 2.500000000000000000e-01
+#Next
+#
+#Validating ch5
+#Covariance = 1.163749077455013324e-04
+#Accumulated difference = 3.385911115917152742e-03
+#Difference in maxima = 1.326867779957524585e-03
+#Difference in minima = 9.334970576141445042e-04
+#Timelag between maxima = 4.999999999999715783e-02
+#Timelag between minima = 7.999999999999989342e-01
+#Next
+#
+#Validating ch7
+#Covariance = 1.126985405540967749e-04
+#Accumulated difference = 2.857119198201266575e-03
+#Difference in maxima = 3.582696745840288632e-04
+#Difference in minima = 4.523664511720201960e-03
+#Timelag between maxima = 0.000000000000000000e+00
+#Timelag between minima = 2.150000000000000355e+00
+#Next
+#
+#Validating ch9
+#Covariance = 1.250280909578863607e-04
+#Accumulated difference = 3.247257369996859894e-03
+#Difference in maxima = 2.212919173501755321e-03
+#Difference in minima = 3.427581218447350343e-03
+#Timelag between maxima = 0.000000000000000000e+00
+#Timelag between minima = 1.150000000000000355e+00
+
+expected_covariance = {'Boundary': 5.269749049737429159e-05, 
+                       'ch5': 1.163749077455013324e-04,
+                       'ch7': 1.126985405540967749e-04,
+                       'ch9': 1.250280909578863607e-04}
+
+expected_difference = {'Boundary': 8.337114775060755880e-04,
+                       'ch5': 3.385911115917152742e-03,
+                       'ch7': 2.857119198201266575e-03,
+                       'ch9': 3.247257369996859894e-03}
+
+expected_maximum_diff = {'Boundary': 7.166795837686146253e-05,
+                         'ch5': 1.326867779957524585e-03,
+                         'ch7': 3.582696745840288632e-04,
+                         'ch9': 2.212919173501755321e-03}
+
+expected_minimum_diff = {'Boundary': 1.335627610594206094e-04,
+                         'ch5': 9.334970576141445042e-04,
+                         'ch7': 4.523664511720201960e-03,
+                         'ch9': 3.427581218447350343e-03}
+
+expected_argmax_timelag = {'Boundary': 3.000000000000007105e-01,
+                           'ch5': 4.999999999999715783e-02,
+                           'ch7': 0.000000000000000000e+00,
+                           'ch9': 0.000000000000000000e+00}
+
+expected_argmin_timelag = {'Boundary': 2.500000000000000000e-01,
+                           'ch5': 7.999999999999989342e-01,
+                           'ch7': 2.150000000000000355e+00,
+                           'ch9': 1.150000000000000355e+00}
+
+
+
+#-------------------------
+# Read validation dataa
+#-------------------------
+
+print 'Reading', project.boundary_filename
+fid = NetCDFFile(project.boundary_filename, 'r')
+input_time = fid.variables['time'][:]
+validation_data['Boundary'] = fid.variables['stage'][:]
+
+reference_time = []
+fid = open(project.validation_filename)
+lines = fid.readlines()
+fid.close()
+
+for i, line in enumerate(lines[1:]):
+    if i == len(input_time): break
+    
+    fields = line.split()
+
+    reference_time.append(float(fields[0]))    # Record reference time
+    for j, key in enumerate(gauge_names[1:]):  # Omit boundary gauge
+        value = float(fields[1:][j])           # Omit time 
+        validation_data[key].append(value/100) # Convert cm2m
+
+
+# Checks
+assert reference_time[0] == 0.0
+assert reference_time[-1] == finaltime
+assert allclose(reference_time, input_time)
+
+for key in gauge_names:
+    validation_data[key] = ensure_numeric(validation_data[key])
+
+
+#--------------------------------------------------
+# Read and interpolate model output
+#--------------------------------------------------
+
+import sys
+if len(sys.argv) > 1:
+    sww_filename = sys.argv[1]
+else:    
+    sww_filename = project.output_filename
+    
+#f = file_function('okushiri_new_limiters.sww',   #The best so far
+#f = file_function('okushiri_as2005_with_mxspd=0.1.sww',
+f = file_function(sww_filename,
+                  quantities='stage',
+                  interpolation_points=gauge_locations,
+                  use_cache=True,
+                  verbose=True)
+
+
+#--------------------------------------------------
+# Compare model output to validation data
+#--------------------------------------------------
+
+eps = get_machine_precision()
+for k, name in enumerate(gauge_names):
+    sqsum = 0
+    denom = 0
+    model = []
+    print 
+    print 'Validating ' + name
+    observed_timeseries = validation_data[name]
+    for i, t in enumerate(reference_time):
+        model.append(f(t, point_id=k)[0])
+
+
+    # Covariance measures    
+    res = cov(observed_timeseries, model)     
+    print 'Covariance = %.18e' %res
+
+    if res < expected_covariance[name]-eps:
+        print '  Result is better than expected by: %.18e'\
+              %(res-expected_covariance[name])
+        print '  Expect = %.18e' %expected_covariance[name]    
+    elif res > expected_covariance[name]+eps:
+        print '  FAIL: Result is worse than expected by: %.18e'\
+                                %(res-expected_covariance[name])
+        print '  Expect = %.18e' %expected_covariance[name]
+    else:
+        pass
+    
+
+    # Difference measures    
+    res = sum(abs(observed_timeseries-model))/len(model)     
+    print 'Accumulated difference = %.18e' %res
+
+    if res < expected_difference[name]-eps:
+        print '  Result is better than expected by: %.18e'\
+              %(res-expected_difference[name])
+        print '  Expect = %.18e' %expected_difference[name]    
+    elif res > expected_difference[name]+eps:
+        print '  FAIL: Result is worse than expected by: %.18e'\
+                                %(res-expected_difference[name])
+        print '  Expect = %.18e' %expected_difference[name]
+    else:
+        pass
+
+
+    # Extrema
+    res = abs(max(observed_timeseries)-max(model))
+    print 'Difference in maxima = %.18e' %res
+
+    if res < expected_maximum_diff[name]-eps:
+        print '  Result is better than expected by: %.18e'\
+              %(res-expected_maximum_diff[name])
+        print '  Expect = %.18e' %expected_maximum_diff[name]    
+    elif res > expected_maximum_diff[name]+eps:
+        print '  FAIL: Result is worse than expected by: %.18e'\
+                                %(res-expected_maximum_diff[name])
+        print '  Expect = %.18e' %expected_maximum_diff[name]
+    else:
+        pass
+    
+
+    
+    res = abs(min(observed_timeseries)-min(model))
+    print 'Difference in minima = %.18e' %res
+
+    if res < expected_minimum_diff[name]-eps:
+        print '  Result is better than expected by: %.18e'\
+              %(res-expected_minimum_diff[name])
+        print '  Expect = %.18e' %expected_minimum_diff[name]    
+    elif res > expected_minimum_diff[name]+eps:
+        print '  FAIL: Result is worse than expected by: %.18e'\
+                                %(res-expected_minimum_diff[name])
+        print '  Expect = %.18e' %expected_minimum_diff[name]
+    else:
+        pass
+    
+
+    # Locations of extrema
+    i0 = argmax(observed_timeseries)
+    i1 = argmax(model)
+    res = abs(reference_time[i1] - reference_time[i0])
+    print 'Timelag between maxima = %.18e' %res
+
+    if res < expected_argmax_timelag[name]-eps:
+        print '  Result is better than expected by: %.18e'\
+              %(res-expected_argmax_timelag[name])
+        print '  Expect = %.18e' %expected_argmax_timelag[name]    
+    elif res > expected_argmax_timelag[name]+eps:
+        print '  FAIL: Result is worse than expected by: %.18e'\
+                                %(res-expected_argmax_timelag[name])
+        print '  Expect = %.18e' %expected_argmax_timelag[name]
+    else:
+        pass
+    
+
+    i0 = argmin(observed_timeseries)
+    i1 = argmin(model)
+    res = abs(reference_time[i1] - reference_time[i0])
+    print 'Timelag between minima = %.18e' %res
+    if res < expected_argmin_timelag[name]-eps:
+        print '  Result is better than expected by: %.18e'\
+              %(res-expected_argmin_timelag[name])
+        print '  Expect = %.18e' %expected_argmin_timelag[name]    
+    elif res > expected_argmin_timelag[name]+eps:
+        print '  FAIL: Result is worse than expected by: %.18e'\
+                                %(res-expected_argmin_timelag[name])
+        print '  Expect = %.18e' %expected_argmin_timelag[name]
+    else:
+        pass
+
+
+
+
+    if plotting is True:
+        ion()
+        hold(False)
+    
+        plot(reference_time, validation_data[name], 'r-',
+             reference_time, model, 'k-')
+        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')
+
+
+
Index: /anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py
===================================================================
--- /anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py	(revision 3915)
+++ /anuga_validation/automated_validation_tests/okushiri_tank_validation/create_okushiri.py	(revision 3916)
@@ -8,4 +8,5 @@
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.coordinate_transforms.geo_reference import Geo_reference
+from anuga.geospatial_data import Geospatial_data
 
 import project
@@ -13,18 +14,22 @@
 
 def prepare_timeboundary(filename):
-    """Converting benchmark 2 time series to NetCDF tms file.
+    """Convert benchmark 2 time series to NetCDF tms file.
     This is a 'throw-away' code taylor made for files like
     'Benchmark_2_input.txt' from the LWRU2 benchmark
     """
 
-    print 'Preparing time boundary from %s' %filename
+    from Scientific.IO.NetCDF import NetCDFFile
     from Numeric import array
 
-    fid = open(filename)
 
-    #Skip first line
+    print 'Creating', filename
+
+    # Read the ascii (.txt) version of this file
+    fid = open(filename[:-4] + '.txt')
+
+    # Skip first line
     line = fid.readline()
 
-    #Read remaining lines
+    # Read remaining lines
     lines = fid.readlines()
     fid.close()
@@ -42,11 +47,7 @@
 
 
-    #Create tms file
-    from Scientific.IO.NetCDF import NetCDFFile
+    # Create tms NetCDF file
 
-    outfile = filename[:-4] + '.tms'
-    print 'Writing to', outfile
-    fid = NetCDFFile(outfile, 'w')
-
+    fid = NetCDFFile(filename, 'w')
     fid.institution = 'Geoscience Australia'
     fid.description = 'Input wave for Benchmark 2'
@@ -71,4 +72,5 @@
 if __name__ == "__main__":
 
+
     # Prepare time boundary
     prepare_timeboundary(project.boundary_filename)
@@ -82,7 +84,8 @@
     #--------------------------------------------------------------------------
 
-    base_resolution = 10
 
-    # Basic geometry
+    base_resolution = 1 # Use this to coarsen or refine entire mesh. 
+
+    # Basic geometry and bounding polygon
     xleft   = 0
     xright  = 5.448
@@ -90,5 +93,4 @@
     ytop    = 3.402
 
-    # Outline 
     point_sw = [xleft, ybottom]
     point_se = [xright, ybottom]
@@ -96,9 +98,11 @@
     point_ne = [xright, ytop]
 
-    # Midway points (left)
-    point_mtop = [xleft + (xright-xleft)/3+1, ytop]
-    point_mbottom = [xleft + (xright-xleft)/3+1, ybottom]
+    bounding_polygon = [point_se,
+                        point_ne,
+                        point_nw,
+                        point_sw]
+    
 
-    #Localised refined area for gulleys
+    # Localised refined area for gulleys
     xl = 4.8
     xr = 5.3
@@ -109,16 +113,19 @@
     p2 = [xr, yt]
     p3 = [xl, yt]
+    
     gulleys = [p0, p1, p2, p3]
+    
 
-    #Island area
+
+    # Island area and drawdown region
     island_0 = [xleft + 2*(xright-xleft)/3+1.2, ytop-0.5]
     island_1 = [xleft + 2*(xright-xleft)/3+0.5, ybottom + 2*(ytop-ybottom)/3]
-    island_2 = [xleft + (xright-xleft)/2+0.3, ybottom + 2*(ytop-ybottom)/3-0.3]
-    island_3 = [xleft + (xright-xleft)/2+0.3, ybottom + (ytop-ybottom)/3+0.3]
+    island_2 = [xleft + (xright-xleft)/2+0.4, ybottom + 2*(ytop-ybottom)/3-0.3]
+    island_3 = [xleft + (xright-xleft)/2+0.4, ybottom + (ytop-ybottom)/3+0.3]
     island_4 = [xleft + 2*(xright-xleft)/3+0.4, ybottom + (ytop-ybottom)/3-0.3]
-    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.2]
-    island_6 = [xl-.01, yb]  #OK
-    island_7 = [xl-.01, yt]  #OK      
-
+    island_5 = [xleft + 2*(xright-xleft)/3+1.2, ybottom+0.8]
+    island_6 = [xl-.01, yb]  # Keep right edge just off the gulleys
+    island_7 = [xl-.01, yt]
+ 
     island = [island_0, island_1, island_2,
               island_3, island_4, island_5,
@@ -126,31 +133,27 @@
 
 
-    
 
-    bounding_polygon = [point_se,
-                        point_ne,
-                        point_mtop,
-                        point_nw,
-                        point_sw,
-                        point_mbottom]
-    
+    # Region spanning half right hand side of domain just inside boundary
+    rhs_nw = [xleft + (xright-xleft)/3+1, ytop-1.4]
+    rhs_sw = [xleft + (xright-xleft)/3+1, ybottom+0.5]
+    rhs_se = [xright-0.1, ybottom+0.2]
+    rhs_ne = [xright-0.1, ytop-0.2]        
 
-    interior_regions = [[gulleys, 0.00002*base_resolution],
-                        [island, 0.00007*base_resolution]]
-    
-
+    rhs_region = [rhs_nw, rhs_ne, rhs_se, rhs_sw]
 
     
-
-    print 'start create mesh from regions'
+    # Interior regions and creation of mesh
+    interior_regions = [[rhs_region, 0.0005],
+                        [island, 0.0003*base_resolution],
+                        [gulleys, 0.00003*base_resolution]]    
 
     meshname = project.mesh_filename + '.msh'
     m = create_mesh_from_regions(bounding_polygon,
-                                 boundary_tags={'wall': [0, 1, 2, 4, 5],
-                                                'wave': [3]},
+                                 boundary_tags={'wall': [0, 1, 3],
+                                                'wave': [2]},     
                                  maximum_triangle_area=0.1*base_resolution,
                                  interior_regions=interior_regions,
                                  filename=project.mesh_filename,
-                                 use_cache=True,
+                                 use_cache=False,
                                  verbose=True)
 
Index: /anuga_validation/automated_validation_tests/okushiri_tank_validation/project.py
===================================================================
--- /anuga_validation/automated_validation_tests/okushiri_tank_validation/project.py	(revision 3915)
+++ /anuga_validation/automated_validation_tests/okushiri_tank_validation/project.py	(revision 3916)
@@ -1,8 +1,21 @@
-"""Common filenames for LWRU benchmark project 2
+"""Common filenames for Okushiri Island validation
+Formats are given as ANUGA native netCDF where applicable.
+
 """
 
-boundary_filename = 'Benchmark_2_input.txt'
-bathymetry_filename = 'Benchmark_2_Bathymetry.txt'
+# Given boundary wave
+boundary_filename = 'Benchmark_2_input.tms'
+
+# Observed timeseries
+validation_filename = 'output_ch5-7-9.txt'
+
+# Digital Elevation Model
+bathymetry_filename = 'Benchmark_2_Bathymetry.pts'
+
+# Triangular mesh
 mesh_filename = 'Benchmark_2.msh'
+
+# Model output
+output_filename = 'okushiri_auta_validation.sww'
 
 
Index: /anuga_validation/automated_validation_tests/okushiri_tank_validation/run_okushiri.py
===================================================================
--- /anuga_validation/automated_validation_tests/okushiri_tank_validation/run_okushiri.py	(revision 3916)
+++ /anuga_validation/automated_validation_tests/okushiri_tank_validation/run_okushiri.py	(revision 3916)
@@ -0,0 +1,93 @@
+"""Validation of the AnuGA implementation of the shallow water wave equation.
+
+This script sets up Okushiri Island benchmark as published at the
+
+THE THIRD INTERNATIONAL WORKSHOP ON LONG-WAVE RUNUP MODELS
+June 17-18 2004
+Wrigley Marine Science Center
+Catalina Island, California
+http://www.cee.cornell.edu/longwave/
+
+
+The validation data was downloaded and made available in this directory
+for convenience but the original data is available at
+http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=2
+where a detailed description of the problem is also available.
+
+
+Run create_okushiri.py to process the boundary condition and build a the
+mesh before running this script.
+
+"""
+
+# Module imports
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Reflective_boundary
+from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
+from anuga.abstract_2d_finite_volumes.util import file_function
+
+from anuga_parallel.parallel_api import myid, numprocs, distribute
+
+import project
+
+
+#-------------------------
+# Create Domain from mesh
+#-------------------------
+domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
+print domain.statistics()
+
+
+#-------------------------
+# Initial Conditions
+#-------------------------
+domain.set_quantity('friction', 0.0)
+domain.set_quantity('stage', 0.0)
+domain.set_quantity('elevation',
+                    filename=project.bathymetry_filename,
+                    alpha=0.02,                    
+                    verbose=True,
+                    use_cache=True)
+
+
+#-------------------------
+# Distribute domain if run in parallel
+#-------------------------
+if numprocs > 1:
+    domain = distribute(domain)
+
+
+#-------------------------
+# Set simulation parameters
+#-------------------------
+domain.set_name(project.output_filename)  # Name of output sww file
+domain.set_default_order(2)               # Apply second order scheme 
+domain.set_all_limiters(0.9)              # Max second order scheme (old lim)
+domain.set_minimum_storable_height(0.001) # Don't store w < 0.001m
+domain.set_maximum_allowed_speed(0.1)     # Allow a little runoff (0.1 is OK)
+
+
+#-------------------------
+# Boundary Conditions
+#-------------------------
+
+# Create boundary function from timeseries provided in file
+function = file_function(project.boundary_filename,
+                         domain, verbose=True)
+
+# Create and assign boundary objects
+Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
+Br = Reflective_boundary(domain)
+domain.set_boundary({'wave': Bts, 'wall': Br})
+
+
+#-------------------------
+# Evolve through time
+#-------------------------
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 0.05, finaltime = 22.5):
+    domain.write_time()
+
+print 'That took %.2f seconds' %(time.time()-t0)
Index: uga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py
===================================================================
--- /anuga_validation/automated_validation_tests/okushiri_tank_validation/validate_okushiri.py	(revision 3915)
+++ 	(revision )
@@ -1,229 +1,0 @@
-"""Automated validation of ANUGA.
-
-This script runs the benchmark on a regular grid and verifies that
-the predicted timeseries at selected gauges are close enought to the observed
-timeseries provided from the experiment.
-
-See anuga_validation/okushiri_2005 for scripts that produce simulations at
-greater resolution.
-
-"""
-
-# Module imports
-from Numeric import array, zeros, Float, allclose
-
-from anuga.utilities.numerical_tools import ensure_numeric
-from anuga.shallow_water import Domain
-from anuga.shallow_water import Reflective_boundary
-from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
-from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
-from anuga.abstract_2d_finite_volumes.util import file_function
-
-import project
-
-try:
-    import pylab
-except:
-    plotting = False
-else:
-    plotting = True
-    
-
-finaltime = 22.5
-timestep = 0.05
-
-#-------------------------
-# Create Domain
-#-------------------------
-#N = 50
-#N = 50
-#points, vertices, boundary = rectangular_cross(N, N/5*3,
-#                                               len1=5.448, len2=3.402)
-#domain = Domain(points, vertices, boundary)
-
-
-print 'Creating domain from', project.mesh_filename
-
-domain = Domain(project.mesh_filename, use_cache=True, verbose=True)
-
-
-domain.set_name('okushiri_automated_validation_varmesh')
-domain.set_default_order(2)
-domain.set_minimum_storable_height(0.001)
-domain.set_maximum_allowed_speed(0) # The default in August 2005
-
-
-# Set old (pre Sep 2006) defaults for limiters
-domain.beta_w      = 0.9
-domain.beta_w_dry  = 0.9
-domain.beta_uh     = 0.9
-domain.beta_uh_dry = 0.9
-domain.beta_vh     = 0.9
-domain.beta_vh_dry = 0.9
-
-domain.check_integrity()
-
-print domain.statistics()
-
-#-------------------------
-# Initial Conditions
-#-------------------------
-domain.set_quantity('friction', 0.0)
-domain.set_quantity('stage', 0.0)
-domain.set_quantity('elevation',
-                    filename=project.bathymetry_filename[:-4] + '.pts',
-                    alpha=0.02,                    
-                    verbose=True,
-                    use_cache=True)
-
-#-------------------------
-# Boundary Conditions
-#-------------------------
-Br = Reflective_boundary(domain)
-
-function = file_function(project.boundary_filename[:-4] + '.tms',
-                         domain,
-                         verbose=True)
-
-Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
-
-domain.set_boundary({'wave': Bts, 'wall': Br})
-
-
-
-
-#-------------------------
-# Set up validation data
-#-------------------------
-gauge_locations = [[0.000, 1.696]]  #Boundary
-gauge_locations += [[4.521, 1.196], [4.521, 1.696], [4.521, 2.196]] #Ch 5-7-9
-
-#Boundary input wave
-filename = project.boundary_filename[:-4] + '.tms'
-print 'Reading', filename
-from Scientific.IO.NetCDF import NetCDFFile
-from Numeric import allclose
-fid = NetCDFFile(filename, 'r')#
-
-input_time = fid.variables['time'][:]
-input_stage = fid.variables['stage'][:]
-
-
-# reference gauge timeseries
-reference_time = []
-ch5 = []
-ch7 = []
-ch9 = []
-fid = open('output_ch5-7-9.txt')
-lines = fid.readlines()
-fid.close()
-#for i, line in enumerate(lines[1:]):
-for line in lines[1:]:
-    fields = line.split()
-
-    reference_time.append(float(fields[0]))
-    ch5.append(float(fields[1])/100)   #cm2m
-    ch7.append(float(fields[2])/100)   #cm2m
-    ch9.append(float(fields[3])/100)   #cm2m
-
-    if fields[0] == str(finaltime): break    
-    
-
-assert reference_time[0] == 0.0
-assert reference_time[-1] == finaltime
-
-reference_gauge_values = [ensure_numeric(input_stage),
-                          ensure_numeric(ch5),
-                          ensure_numeric(ch7),
-                          ensure_numeric(ch9)] 
-
-
-
-predicted_gauge_values = [[], [], [], []]
-
-
-#-------------------------
-# Evolve through time 
-#-------------------------
-import time
-t0 = time.time()
-
-w_max = 0
-for j, t in enumerate(domain.evolve(yieldstep = timestep,
-                                    finaltime = finaltime)):
-    domain.write_time()
-
-    # Record gauge values
-    stage = domain.get_quantity('stage')
-    for i, (x, y) in enumerate(gauge_locations):
-        w = stage.get_values(interpolation_points=[[x,y]])[0]
-        predicted_gauge_values[i].append(w)
-
-        #print 'difference', x, y,\
-        #      predicted_gauge_values[i][j] - reference_gauge_values[i][j]
-
-
-    w = domain.get_maximum_inundation_elevation()
-    x, y = domain.get_maximum_inundation_location()
-    print '  Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y)
-    print
-    
-    if w > w_max:
-        w_max = w
-        x_max = x
-        y_max = y
-
-
-print '**********************************************'
-print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)
-print '**********************************************'
-        
-    
-
-#-------------------------
-# Validate
-#-------------------------
-for i, (x, y) in enumerate(gauge_locations):        
-    predicted_gauge_values[i] = ensure_numeric(predicted_gauge_values[i])
-    predicted_gauge_values[i] = predicted_gauge_values[i][:451]    
-
-    print predicted_gauge_values[i].shape
-    print reference_gauge_values[i].shape
-    #print predicted_gauge_values.shape
-    #print reference_gauge_values.shape    
-    
-    sqsum = sum((predicted_gauge_values[i].flat - \
-                 reference_gauge_values[i].flat)**2)
-    denom = sum(reference_gauge_values[i].flat**2)
-
-    print i, sqsum, sqsum/denom
-
-    msg = 'Rms error is %f' %(sqsum/denom)
-    print msg
-    
-    #assert sqsum/denom < 0.01, msg
-    
-
-    if plotting is True:
-        
-        from pylab import *
-
-        ion()
-        hold(False)
-        plot(reference_time, reference_gauge_values[i], 'r.-',
-             reference_time, predicted_gauge_values[i], 'k-')
-
-        title('Gauge %d (%f,%f)' %(i,x,y))
-        xlabel('time(s)')
-        ylabel('stage (m)')    
-        legend(('Observed', 'Modelled'), shadow=True, loc='upper left')
-        savefig('Gauge_%d.png' %i, dpi = 300)
-        
-        #raw_input('Next')
-        
-
-if plotting is True:
-    show()
-    
-
-print 'That took %.2f seconds' %(time.time()-t0)
