Index: /anuga_work/production/australia_ph2/sydney/run_model.py
===================================================================
--- /anuga_work/production/australia_ph2/sydney/run_model.py	(revision 6287)
+++ /anuga_work/production/australia_ph2/sydney/run_model.py	(revision 6287)
@@ -0,0 +1,142 @@
+"""Run a tsunami inundation scenario for Busselton, WA, Australia.
+
+The scenario is defined by a triangular mesh created from project.polygon, the
+elevation data is compiled into a pts file through build_busselton.py and a
+simulated tsunami is generated through an sts file from build_boundary.py.
+
+Input: sts file (build_boundary.py for respective event)
+       pts file (build_busselton.py)
+       information from project file
+Outputs: sww file stored in project.output_run_time_dir 
+The export_results_all.py and get_timeseries.py is reliant
+on the outputs of this script
+
+Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
+Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+import os
+import time
+
+# Related major packages
+from anuga.interface import create_domain_from_regions
+from anuga.interface import Transmissive_stage_zero_momentum_boundary
+from anuga.interface import Dirichlet_boundary
+from anuga.interface import Reflective_boundary
+from anuga.interface import Field_boundary
+from anuga.interface import create_sts_boundary
+from anuga.interface import csv2building_polygons
+
+from anuga.shallow_water.data_manager import start_screen_catcher
+from anuga.shallow_water.data_manager import copy_code_files
+from anuga.utilities.polygon import read_polygon, Polygon_function
+    
+# Application specific imports
+import project  # Definition of file names and polygons
+
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file. Copy script must be before screen_catcher.
+#------------------------------------------------------------------------------
+copy_code_files(project.output_run, __file__, 
+                os.path.dirname(project.__file__)+os.sep+\
+                project.__name__+'.py' )
+start_screen_catcher(project.output_run, 0, 1)
+
+
+#------------------------------------------------------------------------------
+# 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.
+#------------------------------------------------------------------------------
+print 'Create computational domain'
+
+# Read in boundary from ordered sts file
+event_sts = create_sts_boundary(project.event_sts)
+
+# Reading the landward defined points, this incorporates the original clipping
+# polygon minus the 100m contour.
+landward_boundary = read_polygon(project.landward_boundary)
+
+# Combine sts polyline with landward points
+bounding_polygon_sts = event_sts + landward_boundary
+
+# Number of boundary segments
+N = len(event_sts)-1
+
+# Boundary tags refer to project.landward_boundary
+# 4 points equals 5 segments start at N
+boundary_tags={'back': [N+1,N+2,N+3,N+4, N+5],
+               'side': [N,N+6],
+               'ocean': range(N)}
+
+# Build mesh and domain
+domain = create_domain_from_regions(bounding_polygon_sts,
+                                    boundary_tags=boundary_tags,
+                                    maximum_triangle_area=project.res_poly_all,
+                                    interior_regions=project.interior_regions,
+                                    mesh_filename=project.meshes,
+                                    use_cache=True,
+                                    verbose=True)
+print domain.statistics()
+
+domain.set_name(project.scenario_name)
+domain.set_datadir(project.output_run) 
+domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm
+
+
+#------------------------------------------------------------------------------
+# Setup initial conditions
+#------------------------------------------------------------------------------
+print 'Setup initial conditions'
+
+# Set the initial stage in the offcoast region only
+IC = Polygon_function(project.land_initial_conditions,
+                      default=project.tide,
+                      geo_reference=domain.geo_reference)
+domain.set_quantity('stage', IC, use_cache=True, verbose=True)
+domain.set_quantity('friction', project.friction) 
+domain.set_quantity('elevation', 
+                    filename=project.combined_elevation+'.pts',
+                    use_cache=True,
+                    verbose=True,
+                    alpha=project.alpha)
+
+
+#------------------------------------------------------------------------------
+# Setup boundary conditions 
+#------------------------------------------------------------------------------
+print 'Set boundary - available tags:', domain.get_boundary_tags()
+
+Br = Reflective_boundary(domain)
+Bt = Transmissive_stage_zero_momentum_boundary(domain)
+Bf = Field_boundary(project.event_sts+'.sts',
+                    domain, mean_stage=project.tide,
+                    time_thinning=1,
+                    default_boundary=Bd,
+                    boundary_polygon=bounding_polygon_sts,                    
+                    use_cache=True,
+                    verbose=True)
+
+domain.set_boundary({'back': Br,
+                     'side': Bt,
+                     'ocean': Bf}) 
+
+
+#------------------------------------------------------------------------------
+# Evolve system through time
+#------------------------------------------------------------------------------
+t0 = time.time()
+for t in domain.evolve(yieldstep=project.yieldstep, 
+                       finaltime=project.finaltime,
+                       skip_initial_step=False): 
+    print domain.timestepping_statistics()
+    print domain.boundary_statistics(tags='ocean')
+
+print 'Simulation took %.2f seconds' %(time.time()-t0)
