Index: /anuga_work/production/dampier_2006/project.py
===================================================================
--- /anuga_work/production/dampier_2006/project.py	(revision 3798)
+++ /anuga_work/production/dampier_2006/project.py	(revision 3798)
@@ -0,0 +1,236 @@
+"""Common filenames and locations for topographic data, meshes and outputs.
+"""
+
+from os import sep, environ, getenv, getcwd
+from os.path import expanduser, basename
+#from anuga.utilities.polygon import read_polygon
+import sys
+from anuga.coordinate_transforms.redfearn import\
+     degminsec2decimal_degrees,\
+     convert_from_latlon_to_utm
+
+from time import localtime, strftime
+from anuga.geospatial_data.geospatial_data import *
+
+# Location and naming of scenario data
+state = 'western_australia'
+scenario_dir_name = 'karratha_tsunami_scenario_2005'
+
+basename = 'elevation50m'
+boundary_basename = 'SU-AU'
+
+
+
+# Data output
+codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+
+if sys.platform == 'win32':
+    user = getenv('USERPROFILE')
+else:
+    user = getenv('LOGNAME')
+
+# INUNDATIONHOME is the inundation directory, not the data directory.
+home += sep +'data'
+
+#Derive subdirectories and filenames
+time = strftime('%Y%m%d_%H%M%S',localtime()) #gets time for new dir
+outputtimedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'output'+sep+time+sep
+
+meshdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'meshes'+sep
+datadir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'topographies'+sep
+gaugedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'gauges'+sep
+polygondir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'polygons'+sep
+boundarydir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'boundaries'+sep
+outputdir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'output'+sep
+tidedir = home+sep+state+sep+scenario_dir_name+sep+'anuga'+sep+'tide_data'+sep
+
+#gauge_filename = gaugedir + 'gauge_location_broome.csv'
+
+meshname = meshdir + basename
+outputname = outputtimedir + basename  #Used by post processing
+#boundaryname = outputtimedir + basename  #Used by post processing
+
+
+
+
+# Regions
+
+refzone = 50 
+south = degminsec2decimal_degrees(-20,55,0)
+north = degminsec2decimal_degrees(-20,15,0)
+west = degminsec2decimal_degrees(116,17,0)
+east = degminsec2decimal_degrees(117,10,0)
+
+p0 = [south, degminsec2decimal_degrees(116,32,0)]
+p1 = [south, west]
+p2 = [degminsec2decimal_degrees(-20,23,0), west]
+p3 = [north, degminsec2decimal_degrees(116,45,0)]
+p4 = [north, degminsec2decimal_degrees(117,0,0)]
+p5 = [p2[0], degminsec2decimal_degrees(117,8,0)]
+p6 = [degminsec2decimal_degrees(-20,30,0), east]
+p7 = [degminsec2decimal_degrees(-20,38,0), east]
+p8 = [south, east]
+
+
+bounding_polygon, zone =\
+                  convert_from_latlon_to_utm([p0, p1, p2, p3, p4,
+                                              p5, p6, p7, p8])
+refzone = zone
+
+
+#Interior regions
+
+
+# CIPMA point of interest
+cipma_latitude = -20.588456
+cipma_longitude = 116.771527
+
+
+k0 = [cipma_latitude-0.02, cipma_longitude-0.02]
+k1 = [cipma_latitude-0.02, cipma_longitude+0.02]
+k2 = [cipma_latitude+0.02, cipma_longitude+0.02]
+k3 = [cipma_latitude+0.02, cipma_longitude-0.02]
+
+cipma_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+assert zone == refzone
+
+
+
+#Interior regions
+karratha_south = degminsec2decimal_degrees(-20,44,0)
+karratha_north = degminsec2decimal_degrees(-20,42,0)
+karratha_west = degminsec2decimal_degrees(116,48,0)
+karratha_east = degminsec2decimal_degrees(116,53,30)
+
+k0 = [karratha_south, karratha_west]
+k1 = [karratha_south, karratha_east]
+k2 = [karratha_north, karratha_east]
+k3 = [karratha_north, karratha_west]    
+
+karratha_polygon, zone = convert_from_latlon_to_utm([k0, k1, k2, k3])
+assert zone == refzone
+
+
+#Interior regions
+dampier_south = degminsec2decimal_degrees(-20,40,0)
+dampier_north = degminsec2decimal_degrees(-20,38,10)
+dampier_west = degminsec2decimal_degrees(116,43,0)
+dampier_east = degminsec2decimal_degrees(116,45,0)
+
+d0 = [dampier_south, dampier_west]
+d1 = [dampier_south, dampier_east]
+d2 = [dampier_north, dampier_east]
+d3 = [dampier_north, dampier_west]    
+
+dampier_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Interior regions
+refinery_south = degminsec2decimal_degrees(-20,37,50)
+refinery_north = degminsec2decimal_degrees(-20,36,0)
+refinery_west = degminsec2decimal_degrees(116,44,0)
+refinery_east = degminsec2decimal_degrees(116,46,10)
+
+d0 = [refinery_south, refinery_west]
+d1 = [refinery_south, refinery_east]
+d2 = [refinery_north, refinery_east]
+d3 = [refinery_north, refinery_west]    
+
+refinery_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Interior region around 468899, 7715177:
+#lat (-20, 39, 44.93753), lon (116, 42, 5.09106)
+
+point_south = degminsec2decimal_degrees(-20,39,46)
+point_north = degminsec2decimal_degrees(-20,39,42)
+point_west = degminsec2decimal_degrees(116,42,0)
+point_east = degminsec2decimal_degrees(116,42,10)
+
+d0 = [point_south, point_west]
+d1 = [point_south, point_east]
+d2 = [point_north, point_east]
+d3 = [point_north, point_west]    
+
+point_polygon, zone = convert_from_latlon_to_utm([d0, d1, d2, d3])
+assert zone == refzone
+
+
+#Neils areas around interesting points
+neil1_point1 = [degminsec2decimal_degrees(-20,35,34),
+                degminsec2decimal_degrees(116,45,18)]
+neil1_point2 = [degminsec2decimal_degrees(-20,36,15),
+                degminsec2decimal_degrees(116,46,18)]
+neil1_point3 = [degminsec2decimal_degrees(-20,35,9),
+                degminsec2decimal_degrees(116,47,17)]
+neil1_point4 = [degminsec2decimal_degrees(-20,34,26),
+                degminsec2decimal_degrees(116,46,17)]
+
+neil1_polygon, zone = convert_from_latlon_to_utm([neil1_point1,
+                                                  neil1_point2,
+                                                  neil1_point3,
+                                                  neil1_point4])
+assert zone == refzone
+
+
+
+
+neil2_point1 = [degminsec2decimal_degrees(-20,39,36),
+                degminsec2decimal_degrees(116,41,33)]
+neil2_point2 = [degminsec2decimal_degrees(-20,40,10),
+                degminsec2decimal_degrees(116,42,13)]
+neil2_point3 = [degminsec2decimal_degrees(-20,38,39),
+                degminsec2decimal_degrees(116,43,49)]
+neil2_point4 = [degminsec2decimal_degrees(-20,38,5),
+                degminsec2decimal_degrees(116,43,9)]
+
+neil2_polygon, zone = convert_from_latlon_to_utm([neil2_point1,
+                                                  neil2_point2,
+                                                  neil2_point3,
+                                                  neil2_point4])
+assert zone == refzone
+
+
+
+
+
+#Withnell bay
+wb_point1 = [degminsec2decimal_degrees(-20,35,34),
+                degminsec2decimal_degrees(116,45,18)]
+wb_point2 = [degminsec2decimal_degrees(-20,36,15),
+                degminsec2decimal_degrees(116,46,18)]
+wb_point3 = [degminsec2decimal_degrees(-20,35,9),
+                degminsec2decimal_degrees(116,47,17)]
+wb_point4 = [degminsec2decimal_degrees(-20,34,26),
+                degminsec2decimal_degrees(116,46,17)]
+
+wb_polygon, zone = convert_from_latlon_to_utm([wb_point1, wb_point2,
+                                               wb_point3, wb_point4])
+assert zone == refzone
+
+
+
+
+
+#Larger Withnell bay
+lwb_point1 = [degminsec2decimal_degrees(-20,35,59),
+                degminsec2decimal_degrees(116,42,00)]
+lwb_point2 = [degminsec2decimal_degrees(-20,36,50),
+                degminsec2decimal_degrees(116,46,50)]
+lwb_point3 = [degminsec2decimal_degrees(-20,34,00),
+                degminsec2decimal_degrees(116,47,39)]
+lwb_point4 = [degminsec2decimal_degrees(-20,33,00),
+                degminsec2decimal_degrees(116,42,50)]
+
+lwb_polygon, zone = convert_from_latlon_to_utm([lwb_point1, lwb_point2,
+                                                lwb_point3, lwb_point4])
+                                                      
+assert zone == refzone
+
+
+
+
Index: /anuga_work/production/dampier_2006/run_karratha.py
===================================================================
--- /anuga_work/production/dampier_2006/run_karratha.py	(revision 3798)
+++ /anuga_work/production/dampier_2006/run_karratha.py	(revision 3798)
@@ -0,0 +1,138 @@
+"""Script for running tsunami inundation scenario for Karratha, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.outputtimedir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Reflective_boundary
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+from anuga.geospatial_data.geospatial_data import *
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+# filenames
+meshname = project.meshname+'.msh'
+source_dir = project.boundarydir
+
+
+# creates copy of code in output dir if dir doesn't exist
+if access(project.outputtimedir,F_OK) == 0:
+    mkdir (project.outputtimedir)
+copy (dirname(project.__file__) +sep+ project.__name__+'.py',
+      project.outputtimedir + project.__name__+'.py')
+copy (__file__, project.outputtimedir + basename(__file__))
+print 'project.outputtimedir',project.outputtimedir
+
+
+
+#--------------------------------------------------------------------------
+# 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
+#--------------------------------------------------------------------------
+
+
+interior_regions = [#[project.karratha_polygon, 25000],
+                    #[project.dampier_polygon, 2000],
+                    #[project.refinery_polygon, 2000],
+                    #[project.point_polygon, 2000]]
+                    #[project.cipma_polygon, 20000]]
+                    [project.cipma_polygon, 4000]]     # Caused memory error?    
+                    #[project.neil1_polygon, 4000],
+                    #[project.neil2_polygon, 64000]]
+
+
+
+print 'start create mesh from regions'
+meshname = project.meshname + '.msh'
+create_mesh_from_regions(project.bounding_polygon,
+                         boundary_tags={'back': [7, 8], 'side': [0, 6],
+                                        'ocean': [1, 2, 3, 4, 5]},
+                         maximum_triangle_area=200000,
+                         interior_regions=interior_regions,
+                         filename=meshname,
+                         use_cache=True,
+                         verbose=True)
+
+#-------------------------------------------------------------------------
+# Setup computational domain
+#-------------------------------------------------------------------------
+domain = Domain(meshname, use_cache=True, verbose=True)
+print domain.statistics()
+domain.set_name(project.basename)
+domain.set_datadir(project.outputtimedir)
+
+
+#-------------------------------------------------------------------------
+# Setup initial conditions
+#-------------------------------------------------------------------------
+tide = 0.
+domain.set_quantity('stage', tide)
+domain.set_quantity('friction', 0.0) 
+domain.set_quantity('elevation', 
+                    filename = project.datadir + project.basename + '.pts',
+                    use_cache = False,
+                    verbose = True,
+                    alpha = 0.1)
+
+#-------------------------------------------------------------------------
+# Setup boundary conditions 
+#-------------------------------------------------------------------------
+
+print 'Available boundary tags', domain.get_boundary_tags()
+
+Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
+                   domain, verbose = True)
+Br = Reflective_boundary(domain)
+Bd = Dirichlet_boundary([tide,0,0])
+domain.set_boundary({'back': Br,
+                     'side': Bd,
+                     'ocean': Bf}) 
+
+
+#----------------------------------------------------------------------------
+# Evolve system through time
+#----------------------------------------------------------------------------
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 60, finaltime = 28600):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'ocean')     
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
Index: /anuga_work/production/dampier_2006/run_karratha_parallel.py
===================================================================
--- /anuga_work/production/dampier_2006/run_karratha_parallel.py	(revision 3798)
+++ /anuga_work/production/dampier_2006/run_karratha_parallel.py	(revision 3798)
@@ -0,0 +1,185 @@
+"""Script for running tsunami inundation scenario for Karratha, WA, Australia.
+
+Source data such as elevation and boundary data is assumed to be available in
+directories specified by project.py
+The output sww file is stored in project.outputtimedir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated submarine landslide.
+
+Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
+"""
+#------------------------------------------------------------------------------
+# Import necessary modules
+#------------------------------------------------------------------------------
+
+# Standard modules
+from os import sep
+from os.path import dirname, basename
+from os import mkdir, access, F_OK
+from shutil import copy
+import time
+import sys
+
+
+# Related major packages
+from anuga.shallow_water import Domain
+from anuga.shallow_water import Dirichlet_boundary
+from anuga.shallow_water import File_boundary
+from anuga.shallow_water import Reflective_boundary
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+
+from anuga.geospatial_data.geospatial_data import *
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+from anuga_parallel.parallel_api import distribute, myid, numprocs
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+# filenames
+meshname = project.meshname+'.msh'
+source_dir = project.boundarydir
+
+if myid == 0:
+    # creates copy of code in output dir if dir doesn't exist
+    if access(project.outputtimedir,F_OK) == 0:
+        mkdir (project.outputtimedir)
+    copy (dirname(project.__file__) +sep+ project.__name__+'.py',
+          project.outputtimedir + project.__name__+'.py')
+    copy (__file__, project.outputtimedir + basename(__file__))
+    print 'project.outputtimedir',project.outputtimedir
+
+
+    #--------------------------------------------------------------------------
+    # 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
+    #--------------------------------------------------------------------------
+
+
+    print 'start create mesh from regions'
+    from caching import cache
+    meshname = project.meshname + '_%d.msh' %myid
+    _ = cache(create_mesh_from_regions,
+              project.bounding_polygon,
+              {'boundary_tags': {'back': [7, 8], 'side': [0, 6],
+                                 'ocean': [1, 2, 3, 4, 5]},
+               'maximum_triangle_area': 200000,
+               'filename': meshname},
+              verbose = True,
+              evaluate = False)
+
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    domain = Domain(meshname, use_cache = True, verbose = True)
+    print domain.statistics()
+    domain.set_name(project.basename)
+    domain.set_datadir(project.outputtimedir)
+
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+    tide = 0.
+    domain.set_quantity('stage', tide)
+    domain.set_quantity('friction', 0.0) 
+    domain.set_quantity('elevation', 
+                        filename = project.datadir + project.basename + '.pts',
+                        use_cache = False,
+                        verbose = True,
+                        alpha = 0.1)
+    
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+
+    print 'Available boundary tags', domain.get_boundary_tags()
+    #Bf = File_boundary(source_dir + project.boundary_basename + '.sww', 
+    #                   domain, verbose = True)
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([tide,0,0])
+    domain.set_boundary({'back': Br,
+                         'side': Bd,
+                         'ocean': None}) # Bind this one later
+else:
+    domain = None
+
+
+#--------------------------------------------------------------------------
+# Create the parallel domain
+#--------------------------------------------------------------------------
+if numprocs > 1:
+    domain = distribute(domain, verbose=True)
+
+# Set those boundaries that can't be communicated automatically
+# FIXME: There is a problem with the parallel domain.
+# It seems that boundaries are outside coverage.
+#
+# We get
+#
+#------------------------------------------------
+#Interpolation_function (spatio-temporal) statistics:
+#  Extent:
+#    x in [-60494.651593, 149593.652382], len(x) == 56
+#    y in [-74198.017534, 184623.446077], len(y) == 56
+#    t in [0.000000, 28610.000000], len(t) == 145
+#  Quantities:
+#    stage in [-1.797799, 1.644574]
+#    xmomentum in [-33.219473, 21.571094]
+#    ymomentum in [-32.775090, 73.761256]
+#  Interpolation points (xi, eta): number of points == 411
+#    xi in [-420467.297166, -368252.119852]
+#    eta in [-7652664.149808, -7605199.504709]
+#  Interpolated quantities (over all timesteps):
+#    stage at interpolation points in [inf, inf]
+#    xmomentum at interpolation points in [inf, inf]
+#    ymomentum at interpolation points in [inf, inf]
+#------------------------------------------------
+#
+#Instead of
+#
+#------------------------------------------------
+#Interpolation_function (spatio-temporal) statistics:
+#  Extent:
+#    x in [-60494.651593, 149593.652382], len(x) == 56
+#    y in [-74198.017534, 184623.446077], len(y) == 56
+#    t in [0.000000, 28610.000000], len(t) == 145
+#  Quantities:
+#    stage in [-1.797799, 1.644574]
+#    xmomentum in [-33.219473, 21.571094]
+#    ymomentum in [-32.775090, 73.761256]
+#  Interpolation points (xi, eta): number of points == 503
+#    xi in [4747.280215, 96909.729282]
+#    eta in [9240.613414, 83183.465791]
+#  Interpolated quantities (over all timesteps):
+#    stage at interpolation points in [-0.909515, 1.331521]
+#    xmomentum at interpolation points in [-17.054914, 17.110528]
+#    ymomentum at interpolation points in [-24.691402, 16.625063]
+#------------------------------------------------
+
+
+Bf = File_boundary(source_dir + project.boundary_basename + '.sww',
+                   domain, verbose = True)
+domain.modify_boundary({'ocean': Bf})
+
+
+#----------------------------------------------------------------------------
+# Evolve system through time
+#----------------------------------------------------------------------------
+import time
+t0 = time.time()
+
+for t in domain.evolve(yieldstep = 60, finaltime = 28600):
+    domain.write_time()
+    domain.write_boundary_statistics(tags = 'ocean')     
+    
+print 'That took %.2f seconds' %(time.time()-t0)
+
