Index: production/sydney_2006/run_sydney_smf.py
===================================================================
--- production/sydney_2006/run_sydney_smf.py	(revision 2400)
+++ production/sydney_2006/run_sydney_smf.py	(revision 2407)
@@ -6,31 +6,37 @@
 
 The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data and boundary data obtained from a tsunami simulation done with MOST.
+the elevation data and a simulated submarine landslide.
 
-Ole Nielsen, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
+Ole Nielsen and Duncan Gray, GA - 2005 and Adrian Hitchman and Jane Sexton, GA - 2006
 """
 
-tide = 0       #Australian Height Datum (mean sea level)
 
+#-------------------------------------------------------------------------------# Import necessary modules
+#-------------------------------------------------------------------------------
+
+# Standard modules
 import os
 import time
 
+# Related major packages
+from pyvolution.shallow_water import Domain, Reflective_boundary
+from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from pyvolution.combine_pts import combine_rectangular_points_files 
+from pyvolution.pmesh2domain import pmesh_to_domain_instance
 
-from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
-     Dirichlet_boundary, Time_boundary, Transmissive_boundary
-from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
-     dem2pts
-from pyvolution.pmesh2domain import pmesh_to_domain_instance
-from pyvolution.combine_pts import combine_rectangular_points_files 
-from caching import cache
-import project
-from math import pi, sin
-from smf import slump_tsunami, slide_tsunami, Double_gaussian
-from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 
+# Application specific imports
+import project                 # Definition of file names and polygons
+from smf import slump_tsunami  # Function for submarine mudslide 
 
-# Data preparation
+
+#-------------------------------------------------------------------------------
+# 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
+#-------------------------------------------------------------------------------
+
+# filenames
 coarsedemname = project.coarsedemname
 finedemname = project.finedemname
@@ -38,28 +44,17 @@
 
 # coarse data
-cache(convert_dem_from_ascii2netcdf, coarsedemname, {'verbose': True},
-      dependencies = [coarsedemname + '.asc'],
-      verbose = True)
-      #evaluate = True)
+convert_dem_from_ascii2netcdf(coarsedemname, use_cache=True, verbose=True)
+dem2pts(coarsedemname, use_cache=True, verbose=True)
 
-cache(dem2pts, coarsedemname, {'verbose': True},
-      dependencies = [coarsedemname + '.dem'],      
-      verbose = True)
+# fine data (clipping the points file to smaller area)
+convert_dem_from_ascii2netcdf(finedemname, use_cache=True, verbose=True)
+dem2pts(finedemname,
+        easting_min=project.eastingmin,
+        easting_max=project.eastingmax,
+        northing_min=project.northingmin,
+        northing_max= project.northingmax,
+        use_cache=True, 
+        verbose=True)
 
-# fine data
-cache(convert_dem_from_ascii2netcdf, finedemname, {'verbose': True},
-      dependencies = [finedemname + '.asc'],
-      verbose = True)
-      #evaluate = True)
-
-# clipping the fine data
-cache(dem2pts, finedemname, {'verbose': True,
-      'easting_min': project.eastingmin,
-      'easting_max': project.eastingmax,
-      'northing_min': project.northingmin,
-      'northing_max': project.northingmax},
-      dependencies = [finedemname + '.dem'],
-      #evaluate = True,
-      verbose = True)
 
 # combining the coarse and fine data
@@ -67,37 +62,43 @@
                                  project.coarsedemname + '.pts',
                                  project.combineddemname + '.pts')
-                                 
-# Create Triangular Mesh
-# Overall clipping polygon and interior regions defined in project.py
+
+
+#-------------------------------------------------------------------------------                                 
+# Create the triangular mesh 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
+#-------------------------------------------------------------------------------                                 
+
 from pmesh.create_mesh import create_mesh_from_regions
 
-# for whole region
-interior_res = 5000 # maximal area of per triangle
+interior_res = 5000 
 interior_regions = [[project.harbour_polygon_2, interior_res],
-                    [project.botanybay_polygon_2, interior_res]] 
+                    [project.botanybay_polygon_2, interior_res]]
 
-m = cache(create_mesh_from_regions,
-            project.diffpolygonall,
-            {'boundary_tags': {'bottom': [0], 
-                            'right1': [1], 'right0': [2],
-                            'right2': [3], 'top': [4], 'left1': [5],
-                            'left2': [6], 'left3': [7]},
-            'resolution': 100000,
-            'filename': meshname,           
-            'interior_regions': interior_regions},
-            #evaluate=True,    
-            verbose = True)
 
-# Setup domain
-domain = cache(pmesh_to_domain_instance, (meshname, Domain),
-               dependencies = [meshname],                     
-               verbose = True)
+#FIXME: Fix caching of this one once the mesh_interface is ready
+from caching import cache
+_ = cache(create_mesh_from_regions,
+          project.diffpolygonall,
+          {'boundary_tags': {'bottom': [0], 
+                             'right1': [1], 'right0': [2],
+                             'right2': [3], 'top': [4], 'left1': [5],
+                             'left2': [6], 'left3': [7]},
+           'resolution': 100000,
+           'filename': meshname,           
+           'interior_regions': interior_regions,
+           'UTM': True,
+           'refzone': project.refzone},
+          verbose = True)
 
-# Bring in elevation data
-domain.set_quantity('elevation',
-                     filename = project.combineddemname + '.pts',
-                     use_cache = True,
-                     verbose = True)
-  
+
+#-------------------------------------------------------------------------------                                 
+# Setup computational domain
+#-------------------------------------------------------------------------------                                 
+
+domain = pmesh_to_domain_instance(meshname, Domain,
+                                  use_cache = True,
+                                  verbose = True)
+
 print 'Number of triangles = ', len(domain)
 print 'The extent is ', domain.get_extent()
@@ -105,36 +106,51 @@
 domain.set_name(project.basename)
 domain.set_datadir(project.outputdir)
-domain.store = True
-domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
+domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
 
-# Setup Initial conditions
-t = slump_tsunami(length=30000.0,
-                  depth=400.0,
-                  slope=6.0,
-                  thickness=176.0, 
-                  radius=3330,
-                  dphi=0.23,
-                  x0=project.slump_origin[0], 
-                  y0=project.slump_origin[1], 
-                  alpha=0.0, 
-                  domain=domain)
-domain.set_quantity('stage', t)
+
+#-------------------------------------------------------------------------------
+# Set up scenario (tsunami_source is a callable object used with set_quantity)
+#-------------------------------------------------------------------------------
+
+tsunami_source = slump_tsunami(length=30000.0,
+                               depth=400.0,
+                               slope=6.0,
+                               thickness=176.0, 
+                               radius=3330,
+                               dphi=0.23,
+                               x0=project.slump_origin[0], 
+                               y0=project.slump_origin[1], 
+                               alpha=0.0, 
+                               domain=domain)
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup initial conditions
+#-------------------------------------------------------------------------------
+
+domain.set_quantity('stage', tsunami_source)
 domain.set_quantity('friction', 0)
-   
-# Setup Boundary Conditions
-print domain.get_boundary_tags()
+domain.set_quantity('elevation',
+                    filename = project.combineddemname + '.pts',
+                    use_cache = True,
+                    verbose = True)
+
+
+#-------------------------------------------------------------------------------                                 
+# Setup boundary conditions (all reflective)
+#-------------------------------------------------------------------------------
+
+print 'Available boundary tags', domain.get_boundary_tags()
 
 Br = Reflective_boundary(domain)
-Bt = Transmissive_boundary(domain)
-Bd = Dirichlet_boundary([0,0,0])
-# 10 min square wave starting at 1 min, 6m high
-Bw = Time_boundary(domain=domain,
-                   f=lambda t: [(6<t<606)*6, 0, 0])
-
 domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
                       'right2': Br, 'top': Br, 'left1': Br,
                       'left2': Br, 'left3': Br} )
 
-# Evolve
+
+#-------------------------------------------------------------------------------                                 
+# Evolve system through time
+#-------------------------------------------------------------------------------
+
 import time
 t0 = time.time()
