Index: trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py
===================================================================
--- trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 8881)
+++ trunk/anuga_core/source/anuga/abstract_2d_finite_volumes/generic_boundary_conditions.py	(revision 8882)
@@ -576,5 +576,5 @@
                                domain,
                                quantities=domain.conserved_quantities,
-	                       interpolation_points=self.midpoint_coordinates,
+	                           interpolation_points=self.midpoint_coordinates,
                                time_thinning=time_thinning,
                                time_limit=time_limit,
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/build_elevation.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/build_elevation.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/build_elevation.py	(revision 8882)
@@ -0,0 +1,104 @@
+"""Build the elevation data to run a tsunami inundation scenario 
+for Busselton, WA, Australia.
+
+Input: elevation data from project.py
+Output: pts file stored in project.topographies_dir 
+The run_model.py is reliant on the output of this script.
+
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+from os.path import join
+
+# Related major packages
+from anuga.file_conversion.asc2dem import _convert_dem_from_ascii2netcdf as convert_dem_from_ascii2netcdf
+from anuga.file_conversion.dem2pts import dem2pts
+#from anuga. import start_screen_catcher
+from anuga.utilities.file_utils import copy_code_files
+from anuga.geospatial_data.geospatial_data import Geospatial_data
+from anuga.geometry.polygon import read_polygon
+
+# Application specific imports
+import project
+
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+copy_code_files(project.output_build,__file__, 
+                os.path.dirname(project.__file__)+os.sep+\
+                project.__name__+'.py' )
+#start_screen_catcher(project.output_build)
+
+
+#------------------------------------------------------------------------------
+# Preparation of topographic data
+#
+# Convert ASC 2 DEM 2 PTS using source data and store result in source data
+# Do for coarse and fine data
+# Fine pts file to be clipped to area of interest
+#------------------------------------------------------------------------------
+
+print 'project.bounding_polygon', project.bounding_polygon
+print 'project.combined_elevation_basename', project.combined_elevation_basename
+
+# Create Geospatial data from ASCII files
+geospatial_data = {}
+
+for filename in project.ascii_grid_filenames:
+    absolute_filename = join(project.topographies_folder, filename)
+    convert_dem_from_ascii2netcdf(absolute_filename,
+                                  basename_out=absolute_filename,
+                                  use_cache=True,
+                                  verbose=True)
+    dem2pts(absolute_filename, use_cache=True, verbose=True)
+
+    G_grid = Geospatial_data(file_name=absolute_filename+'.pts',
+                                                verbose=True)
+    print 'Clip geospatial object'
+    geospatial_data[filename] = G_grid.clip(project.bounding_polygon)
+
+# Create Geospatial data from TXT files
+
+for filename in project.point_filenames:
+    absolute_filename = join(project.topographies_folder, filename)
+    G_points = Geospatial_data(file_name=absolute_filename,
+                                                verbose=True)
+    print 'Clip geospatial object'
+    geospatial_data[filename] = G_points.clip(project.bounding_polygon)
+
+#-------------------------------------------------------------------------------
+# Combine, clip and export dataset 
+#-------------------------------------------------------------------------------
+extent_polygons = []
+for extent_polygon_filename in project.extent_polygon_filenames:
+    p = read_polygon(join(project.polygons_folder, extent_polygon_filename))
+    extent_polygons.append(p)
+    
+print 'Add geospatial objects' 
+G = None
+for key in geospatial_data:
+    if key == project.point_filenames[0] or key == project.point_filenames[1]:
+        G += geospatial_data[key]
+    elif key == project.point_filenames[2]:
+        D = geospatial_data[key]
+        D = D.clip_outside(extent_polygons[0])
+        D = D.clip_outside(extent_polygons[1])
+        G += D
+    elif key == project.point_filenames[3]:
+        D = geospatial_data[key]
+        D = D.clip_outside(extent_polygons[2])
+        G += D
+
+print 'Export combined DEM file'
+G.export_points_file(project.combined_elevation + '.pts')
+print 'Do txt version too'
+# Use for comparision in ARC
+G.export_points_file(project.combined_elevation + '.txt')
+
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/build_urs_boundary.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/build_urs_boundary.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/build_urs_boundary.py	(revision 8882)
@@ -0,0 +1,167 @@
+"""
+Replacement function that wraps the legacy function urs2sts()
+and passes in the data from the EventSelection <event>.list file.
+"""
+
+import os
+import os.path
+from time import localtime, strftime, gmtime
+
+import project
+
+from Scientific.IO.NetCDF import NetCDFFile
+import numpy as num
+
+from anuga.file_conversion.urs2sts import urs2sts
+
+
+#-------------------------------------------------------------------------------
+# Get gauges (timeseries of index points)
+#-------------------------------------------------------------------------------
+def get_sts_gauge_data(filename, verbose=False):
+    print 'get_sts_gauge_data: filename=%s' % filename
+    fid = NetCDFFile(filename+'.sts', 'r')      #Open existing file for read
+    permutation = fid.variables['permutation'][:]
+    x = fid.variables['x'][:] + fid.xllcorner   #x-coordinates of vertices
+    y = fid.variables['y'][:] + fid.yllcorner   #y-coordinates of vertices
+    points = num.transpose(num.asarray([x.tolist(), y.tolist()]))
+    time = fid.variables['time'][:] + fid.starttime
+    elevation = fid.variables['elevation'][:]
+       
+    basename = 'sts_gauge'
+    quantity_names = ['stage', 'xmomentum', 'ymomentum']
+    quantities = {}
+    for i, name in enumerate(quantity_names):
+        quantities[name] = fid.variables[name][:]
+
+    #---------------------------------------------------------------------------
+    # Get maximum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+     
+    maxname = 'max_sts_stage.csv'
+    print 'get_sts_gauge_data: maxname=%s' % maxname
+    fid_max = open(os.path.join(project.event_folder, maxname), 'w')
+    fid_max.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_max.write('%d, %.6f, %.6f, %.6f\n' % (index, x[j], y[j], max(stage)))
+      
+    #---------------------------------------------------------------------------
+    # Get minimum wave height throughout timeseries at each index point
+    #---------------------------------------------------------------------------
+
+    minname = 'min_sts_stage.csv'
+    fid_min = open(os.path.join(project.event_folder, minname), 'w')
+    fid_min.write('index, x, y, max_stage \n')    
+    for j in range(len(x)):
+        index = permutation[j]
+        stage = quantities['stage'][:,j]
+        xmomentum = quantities['xmomentum'][:,j]
+        ymomentum = quantities['ymomentum'][:,j]
+
+        fid_min.write('%d, %.6f, %.6f, %.6f\n' %(index, x[j], y[j], min(stage)))
+
+        out_file = os.path.join(project.event_folder,
+                                basename+'_'+str(index)+'.csv')
+        fid_sts = open(out_file, 'w')
+        fid_sts.write('time, stage, xmomentum, ymomentum \n')
+
+        #-----------------------------------------------------------------------
+        # End of the get gauges
+        #-----------------------------------------------------------------------
+        for k in range(len(time)-1):
+            fid_sts.write('%.6f, %.6f, %.6f, %.6f\n'
+                          % (time[k], stage[k], xmomentum[k], ymomentum[k]))
+
+        fid_sts.close()      
+
+    fid.close()
+
+    return quantities,elevation,time
+
+
+##
+# @brief Build boundary STS files from one or more MUX files.
+# @param event_file Name of mux meta-file or single mux stem.
+# @param output_dir Directory to write STS data to.
+# @note 'event_file' is produced by EventSelection.
+def build_urs_boundary(event_file, output_dir):
+    '''Build a boundary STS file from a set of MUX files.'''
+
+    # if we are using an EventSelection multi-mux file
+    if project.multi_mux:
+        # get the mux+weight data from the meta-file (in <boundaries>)
+        mux_event_file = os.path.join(project.event_folder, event_file)
+        print 'using multi-mux file', mux_event_file
+        try:
+            fd = open(mux_event_file, 'r')
+            mux_data = fd.readlines()
+            fd.close()
+        except IOError, e:
+            msg = 'File %s cannot be read: %s' % (mux_event_file, str(e))
+            raise Exception, msg
+        except:
+            raise
+
+        # first line of file is # filenames+weight in rest of file
+        num_lines = int(mux_data[0].strip())
+        mux_data = mux_data[1:]
+        print 'number of sources %d' % num_lines
+
+        # quick sanity check on input mux meta-file
+        if num_lines != len(mux_data):
+            msg = ('Bad file %s: %d data lines, but line 1 count is %d'
+                   % (event_file, len(mux_data), num_lines))
+            raise Exception, msg
+
+        # Create filename and weights lists.
+        # Must chop GRD filename just after '*.grd'.
+        mux_filenames = []
+        for line in mux_data:
+            muxname = line.strip().split()[0]
+            split_index = muxname.index('.grd')
+            muxname = muxname[:split_index+len('.grd')]
+            muxname = os.path.join(project.mux_data_folder, muxname)
+            mux_filenames.append(muxname)
+        
+        mux_weights = [float(line.strip().split()[1]) for line in mux_data]
+
+        # Call legacy function to create STS file.
+        print 'creating sts file'
+        urs2sts(mux_filenames,
+                basename_out=output_dir,
+                ordering_filename=project.urs_order,
+                weights=mux_weights,
+                zone=project.zone,
+                mean_stage=project.tide,
+                verbose=True)
+    else:                           # a single mux stem file, assume 1.0 weight
+        mux_file = os.path.join(project.event_folder, event_file)
+        mux_filenames = [mux_file]
+        print 'using single-mux file', mux_file
+
+        weight_factor = 1.0
+        mux_weights = weight_factor*num.ones(len(mux_filenames), num.float)
+            
+        order_filename = project.urs_order
+
+        print 'reading', order_filename
+        # Create ordered sts file
+        urs2sts(mux_filenames,
+                basename_out=output_dir,
+                ordering_filename=order_filename,
+                weights=mux_weights,
+                mean_stage=project.tide,
+                verbose=True)
+
+    # report on progress so far
+    sts_file = os.path.join(project.event_folder, project.scenario_name)
+    quantities, elevation, time = get_sts_gauge_data(sts_file, verbose=False)
+    print len(elevation), len(quantities['stage'][0,:])
+
+    
+
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/cmpsww.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/cmpsww.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/cmpsww.py	(revision 8882)
@@ -0,0 +1,442 @@
+#!/usr/bin/env python
+
+'''
+A program to compare two SWW  files for "equality".
+
+This program makes lots of assumptions about the structure of the SWW files, 
+so if that structure changes, this program must change.
+'''
+
+import sys
+import os
+import os.path
+import getopt
+from Scientific.IO.NetCDF import NetCDFFile
+import numpy as num
+from anuga.config import netcdf_mode_r
+
+
+#####
+# Various constants.
+#####
+
+# Global for the '-q' quiet flag
+quiet = None
+
+# default tolerances - same as numpy defaults
+default_abs_tolerance = 1.0e-08
+default_rel_tolerance = 1.0000000000000001e-05
+
+# Global attributes that should exist and be same in both files
+# Don't have to have all of these, and we don't care about others.
+expect_global_attributes = ['smoothing', 'vertices_are_stored_uniquely',
+                            'order', 'starttime',
+                            'xllcorner', 'yllcorner',
+                            'zone', 'false_easting', 'false_northing',
+                            'datum', 'projection', 'units']
+
+# dimensions expected, with expected values (None means unknown)
+expected_dimensions = {'number_of_volumes': None,
+                       'number_of_vertices': 3,
+                       'numbers_in_range': 2,
+                       'number_of_points': None,
+                       'number_of_timesteps': None}
+
+# Variables expected, with expected dimensions.
+# Don't have to have all of these, and we don't care about others.
+expected_variables = {'x': ('number_of_points',),
+                      'y': ('number_of_points',),
+                      'elevation': ('number_of_points',),
+                      'elevation_range': ('numbers_in_range',),
+                      'z': ('number_of_points',),
+                      'volumes': ('number_of_volumes', 'number of vertices'),
+                      'time': ('number_of_timesteps',),
+                      'stage': ('number_of_timesteps', 'numbers_of_points',),
+                      'stage_range': ('numbers_in_range',),
+                      'xmomentum': ('number_of_timesteps', 'number_of_points'),
+                      'xmomentum_range': ('numbers_in_range'),
+                      'ymomentum': ('number_of_timesteps', 'number_of_points'),
+                      'ymomentum_range': ('numbers_in_range')}
+
+##
+# @brief An exception to inform user of usage problems.
+class Usage(Exception):
+    def __init__(self, msg):
+        self.msg = msg
+
+
+##
+# @brief Compare two SWW files.
+# @param files A tuple of two filenames.
+# @param globals A list of global attribute names to compare.
+# @param timesteps A list of timesteps to compare at.
+# @param variables A list of variable names to compare.
+# @return Returns if files 'equal', else raises RuntimeError.
+def files_are_the_same(files, globals=None, timesteps=None, variables=None,
+                       rel_tolerance=default_rel_tolerance,
+                       abs_tolerance=default_abs_tolerance):
+    # split out the filenames and check they exist
+    (file1, file2) = files
+    filename1 = os.path.basename(file1)
+    filename2 = os.path.basename(file2)
+#    width = max(len(filename1), len(filename2))
+#    filename1 = filename1.rjust(width)
+#    filename2 = filename2.rjust(width)
+
+    error = False
+    error_msg = ''
+
+    try:
+        fid1 = NetCDFFile(file1, netcdf_mode_r)
+    except:
+        error_msg += "\nFile '%s' can't be opened?\n" % file1
+        error = True
+
+    try:
+        fid2 = NetCDFFile(file2, netcdf_mode_r)
+    except:
+        error_msg += "\nFile '%s' can't be opened?\n" % file2
+        error = True
+        fid1.close()
+
+    if error:
+        raise RuntimeError, error_msg
+
+    if globals is None:
+        globals = expect_global_attributes
+
+    #####
+    # First, check that files have the required structure
+    #####
+
+    # dimensions - only check expected dimensions
+    for key in expected_dimensions:
+        if key not in fid1.dimensions.keys():
+            error_msg += ("\nFile %s doesn't contain dimension '%s'\n"
+                          % (filename1, key))
+            error = True
+        if key not in fid2.dimensions.keys():
+            error_msg += ("\nFile %s doesn't contain dimension '%s'\n"
+                          % (filename2, key))
+            error = True
+
+    # now check that dimensions are the same length
+    # NOTE: DOESN'T CHECK 'UNLIMITED' DIMENSIONS YET! (get None at the moment)
+    for dim in expected_dimensions:
+        dim1_shape = fid1.dimensions.get(dim, None)
+        dim2_shape = fid2.dimensions.get(dim, None)
+        if dim1_shape != dim2_shape:
+            error_msg += ("\nFile %s has '%s' dimension of size %s,\n"
+                          "file %s has that dimension of size %s\n"
+                          % (filename1, dim, str(dim1_shape),
+                             filename2, str(dim2_shape)))
+            error = True
+
+    # check that we have the required globals
+    if globals:
+        for glob in globals:
+            if glob not in dir(fid1):
+                error_msg += ("\nGlobal attribute '%s' isn't in file %s\n"
+                              % (glob, filename1))
+                error = True
+            if glob not in dir(fid2):
+                error_msg += ("\nGlobal attribute '%s' isn't in file %s\n"
+                              % (glob, filename2))
+                error = True
+    else:
+        # get list of global attributes
+        glob_vars1 = []
+        glob_vars2 = []
+        for glob in expect_global_attributes:
+            if glob in dir(fid1):
+                glob_vars1.append(glob)
+            if glob in dir(fid2):
+                glob_vars2.append(glob)
+
+        # now check attribute lists are same
+        if glob_vars1 != glob_vars2:
+            error_msg = ('\nFiles differ in global attributes:\n'
+                         '%s: %s,\n'
+                         '%s: %s\n' 
+                         % (filename1, str(glob_vars1),
+                            filename2, str(glob_vars2)))
+            error = True
+
+    # get variables to test
+    if variables:
+        for var in variables:
+            if var not in fid1.variables.keys():
+                error_msg += ("\nVariable '%s' isn't in file %s\n"
+                              % (var, filename1))
+                error = True
+            if var not in fid2.variables.keys():
+                error_msg += ("\nVariable '%s' isn't in file %s\n"
+                              % (var, filename2))
+                error = True
+    else:
+        # check that variables are as expected in both files
+        var_names1 = []
+        var_names2 = []
+        for var_name in expected_variables:
+            if fid1.variables.has_key(var_name):
+                var_names1.append(var_name)
+            if fid2.variables.has_key(var_name):
+                var_names2.append(var_name)
+    
+        if var_names1 != var_names2:
+            error_msg += ('\nVariables are not the same between files:\n'
+                          '%s variables=%s,\n'
+                          '%s variables=%s\n'
+                          % (filename1, str(var_names1), filename2, str(var_names2)))
+            error = True
+        variables = var_names1
+
+    # get size of time dimension
+    num_timesteps1 = fid1.variables['time'].shape
+    num_timesteps2 = fid2.variables['time'].shape
+    if num_timesteps1 != num_timesteps2:
+        error_msg += ('Files have different number of timesteps:\n'
+                      '%s=%s,\n'
+                      '%s=%s\n'
+                      % (filename1, str(num_timesteps1),
+                         filename2, str(num_timesteps2)))
+        error = True
+
+    num_timesteps = num_timesteps1[0]
+
+    # variable shapes same?
+    for var_name in variables:
+        var1 = fid1.variables[var_name]
+        var2 = fid2.variables[var_name]
+        var1_shape = var1.shape
+        var2_shape = var2.shape
+        if var1_shape != var2_shape:
+            error_msg += ('Files differ in variable %s shape:\n'
+                          '%s: %s,\n'
+                          '%s: %s\n'
+                          % (var_name, filename1, str(var1_shape),
+                                       filename2, str(var2_shape)))
+            error = True
+
+    if error:
+        fid1.close()
+        fid2.close()
+        raise RuntimeError, error_msg
+
+    #####
+    # Now check that actual data values are the same
+    #####
+
+    error_msg = ''
+    glob_vars_bad = {}
+    data_vars_bad = {}
+
+    # check values of global attributes
+    for glob_name in globals:
+        if getattr(fid1, glob_name) != getattr(fid2, glob_name):
+            print("\nFiles differ in global '%s':\n"
+                  "%s: '%s',\n"
+                  "%s: '%s'"
+                  % (glob_name, filename1, str(g1), filename2, str(g2)))
+            glob_vars_bad[glob_name] = glob_vars_bad.get(glob_name, 0) + 1
+            error = True
+
+    # check data variables, be clever with time series data
+    max_rel_difference = -1
+    diff_count = 0
+    for var_name in variables:
+        var_dims = expected_variables[var_name]
+        if (len(var_dims) > 1) and (var_dims[0] == 'number_of_timesteps'):
+            # time series, check by timestep block
+            for t in xrange(num_timesteps):
+                var1 = num.array(fid1.variables[var_name][t,:])
+                var2 = num.array(fid2.variables[var_name][t,:])
+                if not num.allclose(var1, var2,
+                                    rtol=rel_tolerance, atol=abs_tolerance):
+                    error = True
+                    for i in xrange(len(var1)):
+                        if not num.allclose(var1[i], var2[i],
+                                            rtol=rel_tolerance,
+                                            atol=abs_tolerance):
+                            abs_difference = num.abs(var1[i]-var2[i])
+                            max_a_b = num.max(num.abs(var1[i]),
+                                              num.abs(var2[i]))
+                            rel_difference = num.abs(abs_difference/max_a_b)
+
+                            if not quiet:
+                                print('\nFiles differ in variable '
+                                      '%s[%d,%d]:\n'
+                                      '%s: %f\n'
+                                      '%s: %f\n'
+                                      'abs. difference=%e, rel. difference=%e\n'
+                                      % (var_name, t, i,
+                                         filename1, var1[i],
+                                         filename2, var2[i],
+                                         abs_difference,
+                                         rel_difference))
+
+                            if rel_difference > max_rel_difference:
+                                max_rel_difference = rel_difference
+                                max_rel_difference_abs = abs_difference
+                                max_rel_difference_a = var1[i]
+                                max_rel_difference_b = var2[i]
+
+                            data_vars_bad[var_name] = data_vars_bad.get(var_name, 0) + 1
+                            diff_count += 1
+        else:
+            # simple data, check whole thing at once
+            var1 = num.array(fid1.variables[var_name][:])
+            var2 = num.array(fid2.variables[var_name][:])
+            if not num.allclose(var1, var2,
+                                rtol=rel_tolerance, atol=abs_tolerance):
+                for j in xrange(len(var1)):
+                    if not num.allclose(var1[j], var2[j],
+                                          rtol=rel_tolerance, atol=abs_tolerance):
+                        abs_difference = num.abs(var1[j]-var2[j])
+                        max_a_b = num.max(num.abs(var1[j]),
+                                          num.abs(var2[j]))
+                        rel_difference = num.abs(abs_difference/max_a_b)
+
+                        if not quiet:
+                            print('\nFiles differ in variable '
+                                  '%s[%d]:\n'
+                                  '%s: %f\n'
+                                  '%s: %f\n'
+                                  'abs. difference=%e, rel. difference=%e\n'
+                                   % (var_name, j,  
+                                      filename1, var1[j],
+                                      filename2, var2[j],
+                                      abs_difference,
+                                      rel_difference))
+
+                        if rel_difference > max_rel_difference:
+                            max_rel_difference = rel_difference
+                            max_rel_difference_abs = abs_difference
+                            max_rel_difference_a = var1[j]
+                            max_rel_difference_b = var2[j]
+
+                        data_vars_bad[var_name] = data_vars_bad.get(var_name, 0) + 1
+                        diff_count += 1
+                error = True
+
+    #####
+    # close files and signal OK or ERROR
+    #####
+
+    fid1.close()
+    fid2.close()
+
+    if error:
+        error_msg += ('\nNumber of data differences=%d\n'
+                      'Maximum relative data difference=%e\n'
+                      'associated absolute difference=%e\n'
+                      "associated 'a' value=%e\n"
+                      "associated 'b' value=%e\n"
+                      % (diff_count, max_rel_difference, max_rel_difference_abs,
+                         max_rel_difference_a, max_rel_difference_b))
+        error_msg += ('\nglob_vars bad=%s\n' % str(glob_vars_bad))
+        error_msg += ('\ndata_vars bad=%s\n' % str(data_vars_bad))
+        raise RuntimeError, error_msg
+
+    return
+
+
+##
+# @brief Return a usage string.
+def usage():
+    result = []
+    a = result.append
+    a('Usage: %s <options> <file1> <file2>\n' % ProgName)
+    a('where <options> is zero or more of:\n')
+    a('                   -h        print this help\n')
+    a('                   -q        be quiet, print only summary of differences\n')
+    a("                   -a <val>  set absolute threshold of 'equivalent'\n")
+    a("                   -r <val>  set relative threshold of 'equivalent'\n")
+    a('                   -g <arg>  check only global attributes specified\n')
+    a('                             <arg> has the form <globname>[,<globname>[,...]]\n')
+    a('                   -t <arg>  check only timesteps specified\n')
+    a('                             <arg> has the form <starttime>[,<stoptime>[,<step>]]\n')
+    a('                   -v <arg>  check only the named variables\n')
+    a('                             <arg> has the form <varname>[,<varname>[,...]]\n')
+    a('and <file1> and <file2> are two SWW files to compare.\n')
+    a('\n')
+    a('The program exit status is one of:\n')
+    a('   0    the two files are equivalent\n')
+    a('   else the files are not equivalent.')
+    return ''.join(result)
+
+##
+# @brief Print a message to stderr.
+def warn(msg):
+    print >>sys.stderr, msg
+
+
+##
+# @brief 
+# @param argv 
+# @return The status code the program will exit with.
+def main(argv=None):
+    global quiet
+
+    if argv is None:
+        argv = sys.argv
+
+    try:
+        try:
+            opts, args = getopt.getopt(argv[1:], 'hqa:g:r:t:v:',
+                                       ['help', 'globals',
+                                        'variables', 'timesteps'])
+        except getopt.error, msg:
+            raise Usage(msg)
+    except Usage, err:
+        print >>sys.stderr, err.msg
+        print >>sys.stderr, "for help use --help"
+        return 2
+
+    # process options
+    globals = None
+    timesteps = None
+    variables = None
+    quiet = False
+    rel_tolerance = default_rel_tolerance
+    abs_tolerance = default_abs_tolerance
+    for opt, arg in opts:
+        if opt in ('-h', '--help'):
+            print usage()
+            sys.exit(0)
+        elif opt in ('-q', '--quiet'):
+            quiet = True
+        elif opt in ('-a', '--absolute'):
+            abs_tolerance = float(arg)
+        elif opt in ('-r', '--relative'):
+            rel_tolerance = float(arg)
+        elif opt in ('-g', '--globals'):
+            globals = arg.split(',')
+        elif opt in ('-t', '--timesteps'):
+            timesteps = arg.split(',')
+        elif opt in ('-v', '--variables'):
+            variables = arg.split(',')
+
+    # process arguments
+    if len(args) != 2:
+        msg = usage()
+        print 'msg=%s' % msg
+        raise Usage(msg)
+
+    try:
+        files_are_the_same(args, globals=globals,
+                           timesteps=timesteps, variables=variables,
+                           rel_tolerance=rel_tolerance,
+                           abs_tolerance=abs_tolerance)
+    except RuntimeError, msg:
+         print msg
+         return 10
+
+
+if __name__ == "__main__":
+    global ProgName
+
+    ProgName = os.path.basename(sys.argv[0])
+
+    sys.exit(main())
+
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/compare_model_timeseries.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/compare_model_timeseries.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/compare_model_timeseries.py	(revision 8882)
@@ -0,0 +1,246 @@
+"""Compare selected timeseries from model outputs
+
+Usage:
+python compare_model_timeseries.py swwfile1 swwfile2
+
+Return 0 if timeseries are 'close enough', otherwise return non zero error code
+
+Typically swwfile1 would be model output from unit test and swwfile2 would
+be a reference model.
+"""
+
+from anuga.abstract_2d_finite_volumes.util import file_function
+from anuga.utilities.numerical_tools import cov, get_machine_precision
+
+
+import pylab
+import numpy as num
+import sys, os
+
+try:
+    import pylab
+    pylab.hold(False)  # Check if this command can be issued
+except:
+    print 'Could not import pylab'
+    plotting = False
+else:
+    # Create plots as png files
+    plotting = True
+
+
+
+# Model specific data
+gauges = {'g10a': [422233.4, 874380.2],
+          'g10b': [422107.9, 873556.8],
+          'g10c': [421966.8, 872890.3],
+          'g10d': [421606.1, 872106.2],
+          'g11a': [422628.2, 873933.2],
+          #'g11b': [422716.2, 873420.6],
+          'g11c': [422689.1, 872859.8],
+          'g11d': [422408.7, 871940.3],
+          'north': [422572.4, 873992.6],
+          'south': [422004.4, 872014.4]}
+
+
+#---------------------------------------
+
+def usage():
+    print 'Usage:'
+    print 'python compare_model_timeseries.py swwfile1 swwfile2 epsilon'
+
+
+
+def get_timeseries(sww_filename, gauges):
+    """Generate time series for sww file based on gauges
+    """
+
+    gauge_locations = gauges.values()
+    gauge_names = gauges.keys()
+
+    tempfile = 'xyz1234tempfile.sww' # Has to end with sww
+
+    os.system('cp %s %s' % (sww_filename, tempfile))
+    f = file_function(tempfile, 
+                      quantities='stage',
+                      interpolation_points=gauge_locations,
+                      use_cache=True,
+                      verbose=True)
+
+    timevector = f.get_time()
+
+    timeseries = {}
+    for k, name in enumerate(gauge_names):
+        model = timeseries[name] = []
+        
+        for t in timevector:
+            model.append(f(t, point_id=k)[0])
+
+
+    return num.array(timevector), timeseries
+
+
+    
+def compare_timeseries(timevector, 
+                       timeseries1,
+                       timeseries2,
+                       name='',
+                       legend = ('Time series 1', 'Time series 2'),
+                       eps=1.0e-6):
+    """Compare and plot two timeseries
+    """
+
+    
+    timeseries1 = num.array(timeseries1)
+    timeseries2 = num.array(timeseries2) 
+    N = timevector.shape[0]
+    assert timeseries1.shape[0] == N    
+    assert timeseries2.shape[0] == N        
+
+    print 'Testing gauge "%s"' % name
+    print 'epsilon', eps
+    # Covariance measure    
+    res = cov(timeseries1-timeseries2)
+    print '2 norm diff', res        
+    msg = '2-norm of timeseries difference was too large: %e' % res
+    assert res < eps, msg
+     
+    # Maximum norm
+    nominator = max(abs(timeseries1-timeseries2))    
+    denominator = max(abs(timeseries1))
+    if denominator > 0:
+        # Relative measure
+        res = nominator/denominator
+    else:
+        # Absolute measure
+        res = nominator/N
+
+    print 'Max abs diff', res    
+    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
+    assert res < eps, msg
+    
+    
+    nominator = sum(abs(timeseries1-timeseries2))    
+    denominator = sum(abs(timeseries1))
+    if denominator > 0:
+        # Relative measure
+        res = nominator/denominator
+    else:
+        # Absolute measure
+        res = nominator/N
+
+    print 'Sum abs diff', res    
+    msg = '%s: Difference between timeseries was too large: %e' % (name, res)
+    assert res < eps, msg
+
+    # Extrema
+    max1 = max(timeseries1)
+    max2 = max(timeseries2)
+    res = abs(max1-max2)
+    print 'Max diff', res
+    msg = '%s: Difference between maxima was too large: %e' % (name, res)
+    assert res < eps, msg    
+
+    min1 = min(timeseries1)
+    min2 = min(timeseries2)
+    res = abs(min1-min2)    
+    print 'Min diff', res
+    msg = '%s: Difference between minima was too large: %e' % (name, res)
+    assert res < eps, msg
+
+    
+    # Locations of extrema
+    i1 = num.argmax(timeseries1)
+    i2 = num.argmax(timeseries2)
+    
+    res = abs(timevector[i1]-timevector[i2])
+    print 'Max loc diff', res
+    msg = '%s: Difference between location of maxima was too large: %e' % (name, res)
+    assert res < eps, msg
+
+    i1 = num.argmin(timeseries1)
+    i2 = num.argmin(timeseries2)
+    
+    res = abs(timevector[i1]-timevector[i2])    
+    print 'Min loc diff', res    
+    msg = '%s: Difference between location of minima was too large: %e' % (name, res)
+    assert res < eps, msg
+    
+
+    if plotting:
+        # Generate plots
+        #pylab.ion() # No plotting on screen
+        pylab.hold(False)
+    
+        pylab.plot(timevector, timeseries1, 'r-',
+                   timevector, timeseries2, 'k-')
+
+        pylab.title('Gauge %s' % name)
+        pylab.xlabel('time(s)')
+        pylab.ylabel('stage (m)')    
+        pylab.legend(legend, shadow=True, loc='upper left')
+        pylab.savefig(name, dpi = 300)        
+
+        # Error vector
+        #pylab.ion() # No plotting on screen
+        pylab.hold(False)
+        pylab.plot(timevector, timeseries1-timeseries2, 'b-')
+        pylab.title('Gauge %s (difference)' % name)    
+        pylab.xlabel('time(s)')
+        pylab.ylabel('stage difference (m)')  
+        pylab.savefig(name + '_diff', dpi = 300)                  
+
+
+    print 'Gauge "%s" OK' % name
+    print
+    
+
+def compare_all_timeseries(swwfile1, swwfile2, eps=1.0e-6):
+    
+    timevector1, timeseries1 = get_timeseries(swwfile1, gauges)
+    timevector2, timeseries2 = get_timeseries(swwfile2, gauges)    
+
+    msg = 'Time vectors were different in models %s and %s' %\
+          (swwfile1, swwfile2)
+    assert num.allclose(timevector1, timevector2), msg
+    timevector = timevector1
+    
+    # Check that both timeseries exist for all gauges
+    for name in timeseries2:
+        assert name in timeseries1
+    
+    for name in timeseries1:
+        assert name in timeseries2
+            
+    # Compare all timeseries data individually
+    for name in timeseries1:
+        compare_timeseries(timevector, 
+                           timeseries1[name],
+                           timeseries2[name],
+                           name=name,
+                           eps=eps)
+        
+        
+        
+            
+    
+    print 'All gauges OK'
+    
+    
+
+if __name__ == '__main__':
+
+    if len(sys.argv) != 4:
+        usage()
+        sys.exit(-1) 
+
+    eps = float(sys.argv[3]) # Get tolerance to be used in comparisons
+    res = compare_all_timeseries(sys.argv[1], sys.argv[2], eps=eps)
+    
+    #try:
+    #    res = compare_all_timeseries(sys.argv[1], sys.argv[2], eps=eps)
+    #except:
+    #    print 'Failed'
+    #    sys.exit(-1)
+    #else:
+    #    print 'Good', res
+    #    sys.exit(0) 
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/download_data_for_patong_scenario.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/download_data_for_patong_scenario.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/download_data_for_patong_scenario.py	(revision 8882)
@@ -0,0 +1,462 @@
+'''
+Required files are downloaded from the ANUGA servers if they are
+out of date or missing.
+
+Modified from the automated_validation_tests/patong/test_patong_scenarios.py code
+'''
+
+import sys
+import os
+import glob
+import unittest
+import time
+import shutil
+
+from anuga.utilities.system_tools \
+     import get_web_file, untar_file, file_length, get_host_name
+import anuga.utilities.log as log
+
+log.log_filename = './validation.log'
+
+# sourceforge download mirror hosts (must end with '/')
+# try these in turn for each file
+## NOTE: It would be more reliable if we could somehow 'poll' Sourceforge
+##       for a list of mirrors instead of hard-coding a list here.  The only
+##       way to do this at the moment is to 'screen-scrape' the data at
+##       http://apps.sourceforge.net/trac/sourceforge/wiki/Mirrors
+##       but that only gets the main web page for the entity, not the
+##       Sourceforge download mirror server.
+MIRRORS = ['http://transact.dl.sourceforge.net/sourceforge/anuga/',       # au
+           'http://voxel.dl.sourceforge.net/sourceforge/anuga/',          # us
+           'http://superb-west.dl.sourceforge.net/sourceforge/anuga/',    # us
+           'http://jaist.dl.sourceforge.net/sourceforge/anuga/',          # jp
+           'http://dfn.dl.sourceforge.net/sourceforge/anuga/'             # de
+          ]
+
+### for testing
+##MIRRORS = ['http://10.7.64.243/patong_validation_data/']       # local linux box
+
+# URL to hand-get data files, if required
+DATA_FILES_URL = ('http://sourceforge.net/project/showfiles.php?'
+                  'group_id=172848&package_id=319323&release_id=677531')
+
+# sequence of mandatory local data objects
+Mandatory_Data_Objects = ('data.tgz',)
+
+# sequence of optional local data objects.
+# these names must be of the form <scene>.sww.<type>.tgz
+# as code below depends upon it.
+Optional_Data_Objects = (
+                         'patong.sww.TRIAL.tgz',
+                         'patong.sww.BASIC.tgz',
+                         'patong.sww.FINAL.tgz'
+                        )
+#Optional_Data_Objects = ('patong.sww.TRIAL.tgz',)
+
+# Associated tolerances to be used in comparisons (would depend on discretisation errors)
+epsilon = {'patong.sww.TRIAL.tgz': 1.0e-3,
+           'patong.sww.BASIC.tgz': 2.0e-3,
+           'patong.sww.FINAL.tgz': 1.0e-3}
+
+
+# path to the local data directory
+Local_Data_Directory = 'local_data'
+
+# path to the remote data directory
+Remote_Data_Directory = 'remote_data'
+
+# name of stdout catch file for runmodel.py
+RUNMODEL_STDOUT = 'runmodel.stdout'
+
+# text at start of 'output dir' line in RUNMODEL_STDOUT file
+OUTDIR_PREFIX = 'Output directory: '
+
+# Name of SWW file produced by run_model.py
+OUTPUT_SWW = 'patong.sww'
+
+
+def setup():
+    '''Prepare for the validation run.'''
+   
+    pass
+
+
+def refresh_local_data(data_objects, target_dir, mirrors):
+    '''Update local data objects from the server.
+
+    data_objects:   list of files to refresh
+    target_dir:     directory in which to put files
+    mirrors:        list of mirror sites to use
+    
+    Each file has an associated *.digest file used to decide
+    if the local file needs refreshing.
+    
+    Return True if all went well, else False.
+    '''
+
+    # decision function to decide if a file contains HTML
+    def is_html(filename):
+        '''Decide if given file contains HTML.'''
+        
+        fd = open(filename)
+        data = fd.read(1024)
+        fd.close()
+
+        if 'DOCTYPE' in data:
+            return True
+        
+        return False
+
+    
+    # local function to get remote file from one of mirrors
+    def get_remote_from_mirrors(remote, local, auth, mirrors):
+        '''Get 'remote' from one of 'mirrors', put in 'local'.'''
+
+        # Get a unique date+time string to defeat caching.  The idea is to add
+        # this to the end of any URL so proxy sees a different request.
+        cache_defeat = '?' + time.strftime('%Y%m%d%H%M%S')
+
+        # try each mirror when getting file
+        for mirror in mirrors:
+            log.debug('Fetching remote file %s from mirror %s'
+                      % (remote, mirror))
+
+            remote_url = mirror + remote + cache_defeat
+            (result, auth) = get_web_file(remote_url, local, auth=auth)
+            if result and is_html(local)==False:
+                log.debug('Success fetching file %s' % remote)
+                return (True, auth)
+            log.debug('Failure fetching from %s' % mirror)
+            auth = None
+
+        log.debug('Failure fetching file %s' % remote)
+        return (False, auth)            
+                
+
+    # local function to compare contents of two files
+    def files_same(file_a, file_b):
+        '''Compare two files to see if contents are the same.'''
+        
+        fd = open(file_a, 'r')
+        data_a = fd.read()
+        fd.close()
+
+        fd = open(file_b, 'r')
+        data_b = fd.read()
+        fd.close()
+
+        return data_a == data_b
+
+        
+    # local function to update one data object
+    def refresh_object(obj, auth, mirrors):
+        '''Update object 'obj' using authentication tuple 'auth'.
+        
+        Return (True, <updated_auth>) if all went well,
+        else (False, <updated_auth>).
+        '''
+
+        # create local and remote file paths.
+        obj_digest = obj + '.digest'
+        
+        remote_file = os.path.join(Remote_Data_Directory, obj)
+        remote_digest = remote_file + '.digest'
+        
+        local_file = os.path.join(Local_Data_Directory, obj)
+        local_digest = local_file + '.digest'
+        
+        # see if missing either digest or object .tgz
+        if not os.path.exists(local_digest) or not os.path.exists(local_file):
+            # no digest or no object, download both digest and object
+            (res, auth) = get_remote_from_mirrors(obj_digest, local_digest, auth, mirrors)
+            if res:
+                (res, auth) = get_remote_from_mirrors(obj, local_file, auth, mirrors)
+        else:
+            # download object digest to remote data directory
+            (res, auth) = get_remote_from_mirrors(obj_digest, remote_digest, auth, mirrors)
+            if res:
+                if not files_same(local_digest, remote_digest):
+                    # digests differ, refresh object
+                    shutil.move(remote_digest, local_digest)
+                    (res, auth) = get_remote_from_mirrors(obj, local_file, auth, mirrors)
+
+        return (res, auth)
+
+    # create local data directory if required
+    log.debug('Creating local directory: %s' % Local_Data_Directory)
+    if not os.path.exists(Local_Data_Directory):
+        os.mkdir(Local_Data_Directory)
+
+    # clean out remote data copy directory
+    log.debug('Cleaning remote directory: %s' % Remote_Data_Directory)
+    shutil.rmtree(Remote_Data_Directory, ignore_errors=True)
+    os.mkdir(Remote_Data_Directory)
+
+    # success, refresh local files
+    auth = None
+    result = True
+    for data_object in data_objects:
+        log.info("Refreshing file '%s'" % data_object)
+        log.debug('refresh_local_data: getting %s from mirrors, auth=%s'
+                  % (data_object, str(auth)))
+        (res, auth) = refresh_object(data_object, auth, mirrors)
+        log.debug('refresh_local_data: returned (res,auth)=%s,%s'
+                  % (str(res), str(auth)))
+        if res == False:
+            log.info('Refresh of file %s failed.' % data_object)
+            result = False
+            # don't use possibly bad 'auth' again,
+            # some proxies lock out on repeated failures.
+            auth = None
+
+    if result:
+        log.critical('Local data has been refreshed.')
+    else:
+        log.critical('Local data has been refreshed, with one or more errors.')
+    log.critical()
+    return result
+
+
+def can_we_run():
+    '''Decide if we can run with the files we have.
+    
+    Return True if we *can* run, else False.
+
+    Tell user what is happening first, then untar files.
+    '''
+
+    log.critical('Checking if you have the required files to run:')
+
+    # get max width of object name string
+    max_width = 0
+    for obj in Mandatory_Data_Objects:
+        max_width = max(len(obj), max_width)
+    for obj in Optional_Data_Objects:
+        max_width = max(len(obj), max_width)
+
+    # if we don't have *all* mandatory object, can't run
+    have_mandatory_files = True
+    for obj in Mandatory_Data_Objects:
+        obj_path = os.path.join(Local_Data_Directory, obj)
+        if os.path.exists(obj_path):
+            log.info('\t%s  found' % obj.ljust(max_width))
+        else:
+            log.info('\t%s  MISSING AND REQUIRED!' % obj.ljust(max_width))
+            have_mandatory_files = False
+
+    # at least *one* of these must exist
+    have_optional_files = False
+    for obj in Optional_Data_Objects:
+        obj_path = os.path.join(Local_Data_Directory, obj)
+        if os.path.exists(obj_path):
+            have_optional_files = True
+            log.info('\t%s  found' % obj.ljust(max_width))
+        else:
+            log.info('\t%s  MISSING!' % obj.ljust(max_width))
+
+    if not have_mandatory_files or not have_optional_files:
+        log.critical('You must obtain the missing files before you can run '
+                     'this validation.')
+        return False
+
+    log.critical('You have enough required files to run.')
+    log.critical()
+
+    return True
+
+
+def set_environment():
+    # modify environment so we use the local data
+    new_inundationhome = os.path.join(Local_Data_Directory, '')
+    os.environ['INUNDATIONHOME'] = new_inundationhome
+    new_muxhome = os.path.join(Local_Data_Directory, 'data')
+    os.environ['MUXHOME'] = new_muxhome
+
+
+def run_simulation(vtype, sim_obj):
+    '''Run a simulation.
+
+    Returns True if all went well, else False.
+    '''
+    
+    # untar the object
+    tar_path = os.path.join(Local_Data_Directory, sim_obj)
+    log.info('Untarring %s in directory %s ...'
+             % (tar_path, Local_Data_Directory))
+    untar_file(tar_path, target_dir=Local_Data_Directory)
+
+    # modify project.py template
+    log.debug("Creating '%s' version of project.py" % vtype)
+    fd = open('project_template.py', 'r')
+    project = fd.readlines()
+    fd.close()
+
+    new_project = []
+    for line in project:
+        new_project.append(line.replace('#!SETUP!#', vtype.lower()))
+            
+    fd = open('project.py', 'w')
+    fd.write(''.join(new_project))
+    fd.close()
+    
+    # import new project.py
+    import project
+
+    # run the simulation, produce SWW file
+    log.info('Running the simulation ...')
+    cmd = 'python run_model.py > %s' % RUNMODEL_STDOUT
+    log.debug("run_simulation: doing '%s'" % cmd)
+    res = os.system(cmd)
+    log.debug("run_simulation: res=%d" % res)
+
+    # 'unimport' project.py
+    del project
+
+    # check result
+    if res != 0:
+        log.critical('Simulation failed, check log')
+
+    return res == 0
+
+def check_that_output_is_as_expected(expected_sww, valid_sww, epsilon):
+    '''Check that validation output is as required.'''
+
+    # get path to expected SWW file
+    log.critical('Checking that simulation results are as expected ...')
+    local_sww = os.path.join(Local_Data_Directory, valid_sww)
+
+    # get output directory from stdout capture file
+    try:
+        fd = open(RUNMODEL_STDOUT, 'r')
+    except IOError, e:
+        log.critical("Can't open catch file '%s': %s"
+                     % (RUNMODEL_STDOUT, str(e)))
+        return 1
+    lines = fd.readlines()
+    fd.close
+
+    output_directory = None
+    for line in lines:
+        if line.startswith(OUTDIR_PREFIX):
+            output_directory = line.replace(OUTDIR_PREFIX, '', 1)
+            output_directory = output_directory.strip()
+            break
+    if output_directory is None:
+        log.critical("Couldn't find line starting with '%s' in file '%s'"
+                     % (OUTDIR_PREFIX, RUNMODEL_STDOUT))
+        return 1
+
+    log.debug('check_that_output_is_as_expected: output_directory=%s'
+              % output_directory)
+    
+    # compare SWW files here and there
+    new_output_sww = os.path.join(output_directory, expected_sww)
+    #cmd = 'python cmpsww.py %s %s > cmpsww.stdout' % (local_sww, new_output_sww)
+    cmd = 'python compare_model_timeseries.py %s %s %e > compare_model_timeseries.stdout' %\
+          (local_sww, new_output_sww, epsilon)
+    print '-------------------------------------'
+    print cmd
+    print '-------------------------------------'    
+    
+    log.debug("check_that_output_is_as_expected: doing '%s'" % cmd)
+    res = os.system(cmd)
+    log.debug("check_that_output_is_as_expected: res=%d" % res)
+    log.critical()
+    print 'Result', res
+    if res == 0:
+        log.info('Simulation results are as expected.')
+    else:
+        log.critical('Simulation results are NOT as expected.')
+        fd = open('compare_model_timeseries.stdout', 'r')
+        cmp_error = fd.readlines()
+        fd.close()
+        log.critical(''.join(cmp_error))
+
+
+def teardown():
+    '''Clean up after validation run.'''
+
+    log.debug('teardown: called')
+    
+    # remove remote directory and stdout capture file
+    #shutil.rmtree(Remote_Data_Directory, ignore_errors=True)
+    #try:
+    #    os.remove(RUNMODEL_STDOUT)
+    #except OSError:
+    #    pass
+            
+
+################################################################################
+# Mainline - run the simulation, check output.
+################################################################################
+
+# set logging levels
+log.console_logging_level = log.INFO
+log.log_logging_level = log.DEBUG
+
+log.debug('Machine we are running on is "%s"' % get_host_name())
+setup()
+
+# prepare user for what is about to happen
+log.critical('''
+This validation requires a working internet connection to refresh its files.
+You may still run this validation without an internet connection if you have the
+required files.
+
+If you are behind a proxy server you will need to supply your proxy details
+such as the proxy server address and your proxy username and password.  These
+can be defined in one or more of the environment variables:
+    HTTP_PROXY
+    PROXY_USERNAME
+    PROXY_PASSWORD
+if you wish.  If not supplied in environment variables you will be prompted for
+the information.
+''')
+
+
+# make sure local data is up to date
+all_objects = Mandatory_Data_Objects + Optional_Data_Objects
+if not refresh_local_data(all_objects, Local_Data_Directory, MIRRORS):
+    if not can_we_run():
+        log.critical("Can't refresh via the internet and you don't have the "
+                     "required files.")
+        log.critical('Terminating the validation.')
+        log.critical('')
+        log.critical('If you get the missing files from %s' % DATA_FILES_URL)
+        log.critical('then you can try to run the validation again.  Put the '
+                     'files into the directory')
+        log.critical("%s." % Local_Data_Directory)
+        sys.exit(10)
+
+# now untar mandatory objects
+for obj in Mandatory_Data_Objects:
+    tar_path = os.path.join(Local_Data_Directory, obj)
+    log.info('Untarring %s in directory %s ...'
+             % (tar_path, Local_Data_Directory))
+    untar_file(tar_path, target_dir=Local_Data_Directory)
+
+# set required environment variables
+set_environment()
+
+## 6/6/2013 -- Removed this code, as purpose now is just to get data
+## now run what simulations we can and check output is as expected
+
+#for odo in Optional_Data_Objects:
+#    start_time = time.time()
+#
+#    _, vtype, _ = odo.rsplit('.', 2)
+#    vtype = vtype.lower()
+#    log.critical('#' * 72)
+#    log.critical('Running Patong "%s" validation ...' % vtype)
+#    if run_simulation(vtype, odo):
+#        # get SWW names expected and valid, check 'equal'
+#        valid_sww, _ = odo.rsplit('.', 1)
+#        expected_sww, _ = valid_sww.rsplit('.', 1)
+#        check_that_output_is_as_expected(expected_sww, valid_sww, epsilon[odo])
+#
+#    stop_time = time.time()
+#    log.critical('"%s" validation took %.1fs\n\n\n' % (vtype, stop_time - start_time))
+#
+## clean up
+#log.critical('Tearing down ...')
+#teardown()
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/project.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/project.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/project.py	(revision 8882)
@@ -0,0 +1,310 @@
+"""
+This file contains all your file and directory definitions 
+for elevation, meshes and outputs.
+"""
+
+import os
+from anuga.utilities.system_tools import get_user_name, get_host_name
+from time import localtime, strftime, gmtime
+from os.path import join, exists
+from anuga.geometry.polygon import read_polygon
+
+#------------------------------------------------------------------------------
+# Directory setup
+#------------------------------------------------------------------------------
+
+# this section needs to be updated to reflect the modelled community.
+# Note, the user needs to set up the directory system accordingly
+state = 'thailand'
+scenario_name = 'patong'
+scenario_folder = 'patong_tsunami_scenario'
+
+#------------------------------------------------------------------------------
+# Initial Conditions
+#------------------------------------------------------------------------------
+
+# Model specific parameters.
+# One or all can be changed each time the run_model script is executed
+tide = 0.8              # difference between MSL and HAT
+zone = 47               # specify zone of model
+model = 'small'         # different size model 'small' or 'large'
+event_number = 'mux_' + model    # the event number or the mux file name
+alpha = 0.1             # smoothing parameter for mesh
+friction = 0.025         # manning's friction coefficient
+starttime = 0           # start time for simulation
+finaltime = 15000       # final time for simulation
+
+use_buildings = False   # FIXME: 6/6/2013 -- Code needs to be edited to support 'True'
+setup = 'basic'         #'final'         # This can be one of three values
+                        #    trial - coarsest mesh, fast
+                        #    basic - coarse mesh
+                        #    final - fine mesh, slowest
+#------------------------------------------------------------------------------
+# Output filename
+#
+# Your output filename should be unique between different runs on different
+# data.
+# The list of items below will be used to create a file in your output
+# directory.
+# Your user name and time+date will be automatically added.  For example,
+#     [setup, tide, event_number]
+# will result in a filename like
+#     20090212_091046_run_final_0_27283_rwilson
+#------------------------------------------------------------------------------
+
+output_comment = [setup, tide, event_number, 'final']
+
+#------------------------------------------------------------------------------
+# Input Data
+#------------------------------------------------------------------------------
+
+# ELEVATION DATA
+# Used in build_elevation.py
+# Format for ascii grids, as produced in ArcGIS + a projection file
+ascii_grid_filenames = []
+
+# Format for point is x,y,elevation (with header)
+point_filenames = ['patong_10m_grid_msl_Project.txt',    # 10m main grid
+                   'patong_10m_small_grid_for_anuga_sub_Project.txt', # 10m saddle
+                   'patong_bay_1s_grid_Project.txt',   # 1s grid
+                   'andaman_3s_grid_Clip2_Project.txt'] #3s grid
+
+### Add csv header list to all files in point_filenames
+##headerlist = ['x', 'y', 'elevation'] 
+##for f in point_filenames:
+##    add_csv_header(join(topographies_folder, f), headerlist)
+
+extent_polygon_filenames = ['patong_10m.txt',
+                            'saddle_10m.txt',
+                            'patong_1s.txt']
+
+
+# BOUNDING POLYGON - for data clipping and estimate of triangles in mesh
+# Used in build_elevation.py
+# Format for points easting,northing (no header)
+bounding_polygon_filename = 'bounding_polygon_'+ model + '.csv'
+bounding_polygon_maxarea = 100000
+
+# INTERIOR REGIONS -  for designing the mesh
+# Used in run_model.py
+# Format for points easting,northing (no header)                    
+interior_regions_data = [['inundation_area.csv', 75],
+                         ['aoi.csv', 200],
+                         ['aos.csv', 900],
+                         ['sw_coast_original.csv', 7000],
+                         ['building_main.csv', 20],
+                         ['building_main_small.csv', 20],
+                         ['building_main_south.csv', 20],
+                         ['building_saddle.csv', 20]]
+
+
+# LAND - used to set the initial stage/water to be offcoast only
+# Used in run_model.py.  Format for points easting,northing (no header)
+land_initial_conditions_filename = [['initial_conditions_' + model +'.csv', 0]]
+
+# GAUGES - for creating timeseries at a specific point
+# Used in get_timeseries.py.  
+# Format easting,northing,name,elevation (with header)
+gauges_filename = 'gauges.csv'
+
+# BUILDINGS POLYOGN - for elevation of buildings
+# Used in run_model.py
+# Format easting,northing,id,floor (with header)
+building_polygon_filename = 'buildings.csv' 
+
+
+# BUILDINGS EXPOSURE - for identifying inundated houses
+# Used in run_building_inundation.py
+# Format latitude,longitude etc (geographic)
+building_exposure_filename = 'busselton_res_clip.csv' # from NEXIS
+
+# BOUNDING POLYGON - used in build_boundary.py and run_model.py respectively
+# NOTE: when files are put together the points must be in sequence
+# For ease go clockwise!
+# Check the run_model.py for boundary_tags
+
+# Thinned ordering file from Hazard Map (geographic)
+# Format is index,latitude,longitude (with header)
+urs_order_filename = 'urs_order_'+ model +'.csv'
+
+# Landward bounding points
+# Format easting,northing (no header)
+landward_boundary_filename = 'landward_boundary_'+ model +'.csv'
+
+# MUX input filename.
+# If a meta-file from EventSelection is used, set 'multi-mux' to True.
+# If a single MUX stem filename (*.grd) is used, set 'multi-mux' to False.
+mux_input_filename = event_number # to be found in event_folder
+                                    # (ie boundaries/event_number/)
+multi_mux = False
+##mux_input_filename = 'event.list'
+##multi_mux = True
+
+#------------------------------------------------------------------------------
+# Clipping regions for export to asc and regions for clipping data
+# Final inundation maps should only be created in regions of the finest mesh
+#------------------------------------------------------------------------------
+
+#CBD extract ascii grid - coordinates from patong_1s extent
+xminCBD = 417445.1119
+xmaxCBD = 425601.7881
+yminCBD = 870663.4547
+ymaxCBD = 876965.3856
+
+
+
+################################################################################
+################################################################################
+####         NOTE: NOTHING WOULD NORMALLY CHANGE BELOW THIS POINT.          ####
+################################################################################
+################################################################################
+
+# Get system user and host names.
+# These values can be used to distinguish between two similar runs by two
+# different users or runs by the same user on two different machines.
+user = get_user_name()
+host = get_host_name()
+
+#-------------------------------------------------------------------------------
+# Output Elevation Data
+#-------------------------------------------------------------------------------
+
+# Output filename for elevation
+# this is a combination of all the data generated in build_elevation.py
+combined_elevation_basename = scenario_name + '_combined_elevation_' + model
+
+#-------------------------------------------------------------------------------
+# Directory Structure
+#-------------------------------------------------------------------------------
+
+# determines time for setting up output directories
+time = strftime('%Y%m%d_%H%M%S', localtime()) 
+gtime = strftime('%Y%m%d_%H%M%S', gmtime()) 
+build_time = time + '_build_' + model
+run_time = time + '_run_'
+
+# create paths for data files.
+output_dirname = os.path.dirname(__file__)
+home = os.path.join(output_dirname, 'local_data', 'data')
+muxhome = home # FIXME (Ole): Get rid off
+
+    
+# check various directories/files that must exist
+anuga_folder = join(home, state, scenario_folder, 'anuga')
+topographies_folder = join(anuga_folder, 'topographies')
+polygons_folder = join(anuga_folder, 'polygons')
+boundaries_folder = join(anuga_folder, 'boundaries')
+output_folder = join(anuga_folder, 'outputs')
+gauges_folder = join(anuga_folder, 'gauges')
+meshes_folder = join(anuga_folder, 'meshes')
+event_folder = join(boundaries_folder, str(event_number))
+
+# MUX data files
+# Directory containing the MUX data files to be used with EventSelection.
+mux_data_folder = join(muxhome, 'mux')
+
+#------------------------------------------------------------------------------
+# Location of input and output data
+#------------------------------------------------------------------------------
+
+# Convert the user output_comment to a string for run_model.py
+output_comment = ('_'.join([str(x) for x in output_comment if x != user])
+                  + '_' + user)
+
+# The absolute pathname of the all elevation, generated in build_elevation.py
+combined_elevation = join(topographies_folder, combined_elevation_basename)
+
+# The absolute pathname of the mesh, generated in run_model.py
+meshes = join(meshes_folder, scenario_name) + model +'.msh'
+
+# The pathname for the urs order points, used within build_urs_boundary.py
+urs_order = join(boundaries_folder, urs_order_filename)
+
+# The absolute pathname for the landward points of the bounding polygon,
+# Used within run_model.py)
+landward_boundary = join(boundaries_folder, landward_boundary_filename)
+
+# The absolute pathname for the .sts file, generated in build_boundary.py
+event_sts = join(event_folder, scenario_name)
+
+# The absolute pathname for the output folder names
+# Used for build_elevation.py
+output_build = join(output_folder, build_time) + '_' + str(user) 
+# Used for run_model.py
+output_run = join(output_folder, run_time) + output_comment 
+# Used by post processing
+output_run_time = join(output_run, scenario_name) 
+
+# The absolute pathname for the gauges file
+# Used for get_timeseries.py
+gauges = join(gauges_folder, gauges_filename)
+
+# The absolute pathname for the gauges file
+# Used for run_model.py
+building_polygon = join(polygons_folder, building_polygon_filename)
+
+# The absolute pathname for the building file
+# Used for run_building_inundation.py
+building_exposure = join(gauges_folder, building_exposure_filename)
+
+# full path to where MUX files (or meta-files) live
+mux_input = join(event_folder, mux_input_filename)
+
+###########################################################################
+#
+# Code below here used to be in setup_model.py
+#
+###########################################################################
+
+sanity_error=False
+
+if setup == 'trial':
+    print 'trial'
+    scale_factor = 100
+    yieldstep = 30.
+elif setup == 'basic': 
+    print 'basic'
+    scale_factor = 4
+    yieldstep = 30.
+elif setup == 'final': 
+    print 'final'
+    scale_factor = 1
+    yieldstep = 30.
+else:
+    print ("Sorry, you must set the 'setup' variable to one of:"
+           '   trial - coarsest mesh, fast\n'
+           '   basic - coarse mesh\n'
+           '   final - fine mesh, slowest\n'
+           '\n'
+           "'setup' was set to '%s'" % setup)
+    sanity_error = True
+
+#-------------------------------------------------------------------------------
+# Check for errors detected above.
+#-------------------------------------------------------------------------------
+
+if sanity_error:
+    msg = 'You must fix the above errors before continuing.'
+    raise Exception, msg
+
+#-------------------------------------------------------------------------------
+# Reading polygons and creating interior regions
+#-------------------------------------------------------------------------------
+
+# Create list of land polygons with initial conditions
+land_initial_conditions = []
+for filename, MSL in land_initial_conditions_filename:
+    polygon = read_polygon(join(polygons_folder, filename))
+    land_initial_conditions.append([polygon, MSL])
+
+# Create list of interior polygons with scaling factor
+interior_regions = []
+for filename, maxarea in interior_regions_data:
+    polygon = read_polygon(join(polygons_folder, filename))
+    interior_regions.append([polygon,
+                                     maxarea*scale_factor])
+
+# Initial bounding polygon for data clipping 
+bounding_polygon = read_polygon(join(polygons_folder,
+                                             bounding_polygon_filename))
+bounding_maxarea = bounding_polygon_maxarea*scale_factor
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/run_model.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/run_model.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/run_model.py	(revision 8882)
@@ -0,0 +1,258 @@
+"""Run a tsunami inundation scenario for Patong, as in Jakeman et al (2010)
+
+The scenario is defined by a triangular mesh created from project.polygon, the
+elevation data is compiled into a pts file through build_elevation.py and a
+simulated tsunami is generated through an sts file from build_boundary.py.
+
+Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
+Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
+Gareth Davies -- 2013
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import os.path
+import time
+from time import localtime, strftime, gmtime
+
+# Related major packages
+from Scientific.IO.NetCDF import NetCDFFile
+
+import anuga
+
+import project
+import anuga.utilities.log as log
+
+# Parallel routines
+from anuga_parallel import distribute, myid, numprocs, finalize, barrier
+
+# Application specific imports
+import build_urs_boundary as bub
+
+#-------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file. Copy script must be before screen_catcher
+#-------------------------------------------------------------------------------
+if(myid==0):
+    # Make output dir and set log filename before anything is logged
+    os.mkdir(project.output_run)
+
+# Tell log module to store log file in output dir
+log.log_filename = os.path.join(project.output_run, 'anuga.log')
+
+if(myid==0):
+    log.critical('Log filename: %s' % log.log_filename)
+    
+    
+    output_dirname = os.path.dirname(project.__file__)
+    anuga.copy_code_files(project.output_run,
+                    [__file__, 
+                     os.path.join(output_dirname, project.__name__+'.py')],
+                    verbose=True
+                    )
+
+
+    #-------------------------------------------------------------------------------
+    # Create the computational domain based on overall clipping polygon with
+    # a tagged boundary and interior regions defined in project.py along with
+    # resolutions (maximal area of per triangle) for each polygon
+    #-------------------------------------------------------------------------------
+
+    log.critical('Create computational domain')
+
+    # Create the STS file
+    # FIXME (Ole): This is deadly dangerous if buildcode changes (as was the case 24th March 2009)
+    # We need to use caching instead!
+    log.critical( 'project.mux_data_folder=%s' % project.mux_data_folder)
+    if not os.path.exists(project.event_sts + '.sts'):
+        bub.build_urs_boundary(project.mux_input_filename, project.event_sts)
+else:
+    print 'Hello from processor ', myid
+
+barrier()
+    
+
+# Read in boundary from ordered sts file
+event_sts = anuga.create_sts_boundary(project.event_sts)
+
+# Reading the landward defined points, this incorporates the original clipping
+# polygon minus the 100m contour
+landward_boundary = anuga.read_polygon(project.landward_boundary)
+
+# Combine sts polyline with landward points
+bounding_polygon_sts = event_sts + landward_boundary
+
+# Number of boundary segments
+num_ocean_segments = len(event_sts) - 1
+# Number of landward_boundary points
+num_land_points = anuga.file_length(project.landward_boundary)
+
+if(myid==0):
+    # Boundary tags refer to project.landward_boundary
+    # 4 points equals 5 segments start at N
+    boundary_tags={'back': range(num_ocean_segments+1,
+                             num_ocean_segments+num_land_points),
+               'side': [num_ocean_segments,
+                        num_ocean_segments+num_land_points],
+               'ocean': range(num_ocean_segments)}
+
+    # Build mesh and domain
+    domain = anuga.create_domain_from_regions(bounding_polygon_sts,
+                                    boundary_tags=boundary_tags,
+                                    maximum_triangle_area=project.bounding_maxarea,
+                                    interior_regions=project.interior_regions,
+                                    mesh_filename=project.meshes,
+                                    use_cache=False,
+                                    verbose=True)
+    log.critical(domain.statistics())
+
+    # FIXME(Ole): How can we make this more automatic?
+    domain.geo_reference.zone = project.zone
+
+
+    domain.set_name(project.scenario_name)
+    domain.set_datadir(project.output_run) 
+    domain.set_flow_algorithm('tsunami')
+    domain.minimum_allowed_height=0.001
+
+    #-------------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------------
+
+    log.critical('Setup initial conditions')
+
+    # Set the initial stage in the offcoast region only
+    if project.land_initial_conditions:
+        IC = anuga.Polygon_function(project.land_initial_conditions,
+                              default=project.tide,
+                              geo_reference=domain.geo_reference)
+    else:
+        IC = project.tide
+
+    domain.set_quantity('stage', IC, use_cache=False, verbose=True)
+    domain.set_quantity('friction', project.friction) 
+    domain.set_quantity('elevation', 
+                        filename=project.combined_elevation+'.pts',
+                        use_cache=False,
+                        verbose=True,
+                        alpha=project.alpha)
+
+    if project.use_buildings:
+        # Add buildings from file
+        log.critical('Reading building polygons')
+        
+        err_mess='ERROR: csv2building_polygons no longer exists (6/6/2013).\
+                  This needs to be fixed before buildings can be ran'
+        raise Exception(err_mess)
+        
+        
+        #building_polygons, building_heights = csv2building_polygons(project.building_polygon)
+        #clipping_polygons=project.building_area_polygons)
+
+        log.critical('Creating %d building polygons' % len(building_polygons))
+        def create_polygon_function(building_polygons, geo_reference=None):
+            L = []
+            for i, key in enumerate(building_polygons):
+                if i%100==0: log.critical(i)
+                poly = building_polygons[key]
+                elev = building_heights[key]
+                L.append((poly, elev))
+                
+                buildings = anuga.Polygon_function(L, default=0.0,
+                                             geo_reference=geo_reference)
+            return buildings
+
+        log.critical('Creating %d building polygons' % len(building_polygons))
+        buildings = cache(create_polygon_function,
+                          building_polygons,
+                          {'geo_reference': domain.geo_reference},
+                          verbose=True)
+
+        log.critical('Adding buildings')
+        domain.add_quantity('elevation',
+                            buildings,
+                            use_cache=False,
+                            verbose=True)
+
+
+else:
+    domain=None
+    print 'Hello from Processor ', myid
+
+barrier()
+
+
+domain=distribute(domain) 
+   
+#-------------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------------
+
+log.critical('Set boundary - available tags:' % domain.get_boundary_tags())
+
+Br = anuga.Reflective_boundary(domain)
+Bs = anuga.Transmissive_stage_zero_momentum_boundary(domain)
+Bf = anuga.Field_boundary(project.event_sts+'.sts',
+                    domain,
+                    mean_stage=project.tide,
+                    time_thinning=20,
+                    default_boundary=anuga.Dirichlet_boundary([0, 0, 0]),
+                    boundary_polygon=bounding_polygon_sts,                    
+                    use_cache=False,
+                    verbose=True)
+
+domain.set_boundary({'back': Br,
+                     'side': Bs,
+                     'ocean': Bf}) 
+
+#-------------------------------------------------------------------------------
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
+t0 = time.time()
+
+## FIXME: GD June 6 2013 -- I believe ANUGA can no longer do multiple evolves
+## as the 'time' variable gets confused. Think this is because 'starttime'
+## exists in both the sts input files, and in the model, and the two are
+## conflated in the code
+##
+## Skip over the first 6000 seconds
+#for t in domain.evolve(yieldstep=2000,
+#                       finaltime=6000):
+#    log.critical(domain.timestepping_statistics())
+#    log.critical(domain.boundary_statistics(tags='ocean'))
+#
+## Start detailed model
+#for t in domain.evolve(yieldstep=project.yieldstep,
+#                       finaltime=project.finaltime,
+#                       skip_initial_step=True):
+#    log.critical(domain.timestepping_statistics())
+#    log.critical(domain.boundary_statistics(tags='ocean'))
+
+
+barrier()
+## Check various important parameters
+import time
+time.sleep(myid)
+print 'Processor ', myid , ' : ', project.yieldstep, domain.flow_algorithm 
+barrier()
+print 'project.yieldstep', project.yieldstep
+for t in domain.evolve(yieldstep=project.yieldstep,
+                       finaltime=project.finaltime):
+    time.sleep(myid*0.01) # To get printing separated
+    print 'Processor ', myid
+    log.critical(domain.timestepping_statistics())
+    log.critical(domain.boundary_statistics(tags='ocean'))
+    print 'Processor ', myid, ' project.yieldstep ', project.yieldstep
+
+## Final print out
+barrier()
+if(myid==0):    
+     log.critical('Simulation took %.2f seconds' %(time.time()-t0))
+
+
+      
Index: trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/validate_patong_scenario.py
===================================================================
--- trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/validate_patong_scenario.py	(revision 8882)
+++ trunk/anuga_core/source/anuga_validation_tests/Case_studies/Patong/patong_run_code/validate_patong_scenario.py	(revision 8882)
@@ -0,0 +1,71 @@
+"""Automatic verification of ANUGA using Patong Beach validation scenario.
+See functions exercised by this wrapper for more details and also the publication (TBA).
+"""
+
+import unittest
+import os
+from subprocess import Popen, PIPE
+
+
+def get_free_memory():
+    """Get available memory. Linux only.
+    """
+    
+    p = Popen('free', shell=True,
+              stdin=PIPE, stdout=PIPE, stderr=PIPE, close_fds=True)
+              
+    if p.stdout is not None:
+        for line in p.stdout.readlines():
+            if line.startswith('Mem'):
+                fields = line.split()
+                mem = int(fields[1]) # Total memory
+    else:
+        mem = None
+        
+        
+    return mem
+
+
+class Test_flow(unittest.TestCase):
+    def setUp(self):
+        for file in os.listdir('.'):    
+            if file.endswith('.stdout') or\
+                    file.endswith('.sww') or\
+                    file.endswith('.msh'):
+                os.remove(file)
+                
+        
+    def tearDown(self):
+        pass
+
+    def test_patong_validation(self):
+        """Exercise Patong Validation Scenario for three resolutions and
+        compare timeseries from modelled results against reference model.
+        """
+        
+        # Bail out if there is insufficient memory.
+        # This works only on *nix systems
+        mem = get_free_memory()
+        if mem is None:
+            msg = 'No information about available memory: '
+            msg += 'Skipping Patong beach validation'
+            raise Exception(msg)
+            
+        if mem < 2000000:
+            msg = 'Insufficient memory to run Patong beach validation. '
+            msg += 'Got %i need at least 8GB. ' % mem
+            msg += 'Skipping Patong beach validation'
+            raise Exception(msg)        
+        
+        #print
+        s = 'test_patong_scenario.py'
+        #print s
+        res = os.system('python %s > test_patong_scenario.stdout' % s)
+        #assert res == 0
+
+
+#-------------------------------------------------------------
+if __name__ == '__main__':
+    suite = unittest.makeSuite(Test_flow, 'test')
+    runner = unittest.TextTestRunner(verbosity=2)
+    runner.run(suite)
