Index: /anuga_work/production/geraldton/build_geraldton.py
===================================================================
--- /anuga_work/production/geraldton/build_geraldton.py	(revision 5007)
+++ /anuga_work/production/geraldton/build_geraldton.py	(revision 5007)
@@ -0,0 +1,163 @@
+"""Script for running tsunami inundation scenario for Dampier, 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.output_time_dir 
+
+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.data_manager import convert_dem_from_ascii2netcdf, dem2pts
+from anuga.geospatial_data.geospatial_data import *
+from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+
+# Application specific imports
+import project   # Definition of file names and polygons
+
+#------------------------------------------------------------------------------
+# Copy scripts to time stamped output directory and capture screen
+# output to file
+#------------------------------------------------------------------------------
+
+start_screen_catcher(project.output_build_time_dir)
+
+print 'time stamp: ',project.gtime
+print 'USER: ', project.user
+
+
+copy_code_files(project.output_build_time_dir,__file__, 
+               dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+
+
+#-------------------------------------------------------------------------------
+# 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
+#-------------------------------------------------------------------------------
+
+# topography directory filenames
+onshore_in_dir_name = project.onshore_in_dir_name
+coast_in_dir_name = project.coast_in_dir_name
+#island_in_dir_name = project.island_in_dir_name
+offshore_in_dir_name = project.offshore_in_dir_name
+
+onshore_dir_name = project.onshore_dir_name
+coast_dir_name = project.coast_dir_name
+#island_dir_name = project.island_dir_name
+offshore_dir_name = project.offshore_dir_name
+
+# creates DEM from asc data
+print "creates DEMs from ascii data"
+convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
+#convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
+
+#creates pts file for onshore DEM
+print "creates pts file for onshore DEM"
+dem2pts(onshore_dir_name,
+#        easting_min=project.eastingmin,
+#        easting_max=project.eastingmax,
+#        northing_min=project.northingmin,
+#        northing_max= project.northingmax,
+        use_cache=True, 
+        verbose=True)
+
+#creates pts file for island DEM
+#dem2pts(island_dir_name, use_cache=True, verbose=True)
+
+print'create Geospatial data1 objects from topographies'
+G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
+G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
+#G3 = Geospatial_data(file_name = island_dir_name + '.pts')
+G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
+
+print'add all geospatial objects'
+G = G1 + G2 + G_off 
+
+print'clip combined geospatial object by bounding polygon'
+#G_clipped = G.clip(project.bounding_polygon)
+#FIXME: add a clip function to pts
+#print'shape of clipped data', G_clipped.get_data_points().shape
+
+print'export combined DEM file'
+if access(project.topographies_dir,F_OK) == 0:
+    mkdir (project.topographies_dir)
+G.export_points_file(project.combined_dir_name + '.xya')
+#G_clipped.export_points_file(project.combined_dir_name + '.xya')
+
+'''
+print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
+G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
+print'split'
+G_all_1, G_all_2 = G_all.split(.10)
+print'export 1'
+G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
+print'export 2'
+G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
+
+
+#-------------------------------------------------------------------------
+# Convert URS to SWW file for boundary conditions 
+#-------------------------------------------------------------------------
+print 'starting to create boundary conditions'
+boundaries_in_dir_name = project.boundaries_in_dir_name
+
+from anuga.shallow_water.data_manager import urs2sww
+
+print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
+print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
+
+#import sys; sys.exit()
+
+#if access(project.boundaries_dir,F_OK) == 0:
+#    mkdir (project.boundaries_dir)
+
+from caching import cache
+cache(urs2sww,
+      (boundaries_in_dir_name,
+       project.boundaries_dir_name1), 
+      {'verbose': True,
+       'minlat': project.south_boundary,
+       'maxlat': project.north_boundary,
+       'minlon': project.west_boundary,
+       'maxlon': project.east_boundary,
+#       'minlat': project.south,
+#       'maxlat': project.north,
+#       'minlon': project.west,
+#       'maxlon': project.east,
+       'mint': 0, 'maxt': 40000,
+#       'origin': domain.geo_reference.get_origin(),
+       'mean_stage': project.tide,
+#       'zscale': 1,                 #Enhance tsunami
+       'fail_on_NaN': False},
+       verbose = True,
+       )
+#       dependencies = source_dir + project.boundary_basename + '.sww')
+
+'''
+
+
+
+
+
+
+
Index: uga_work/production/geraldton/build_gerladton.py
===================================================================
--- /anuga_work/production/geraldton/build_gerladton.py	(revision 5006)
+++ 	(revision )
@@ -1,163 +1,0 @@
-"""Script for running tsunami inundation scenario for Dampier, 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.output_time_dir 
-
-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.data_manager import convert_dem_from_ascii2netcdf, dem2pts
-from anuga.geospatial_data.geospatial_data import *
-from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
-
-# Application specific imports
-import project   # Definition of file names and polygons
-
-#------------------------------------------------------------------------------
-# Copy scripts to time stamped output directory and capture screen
-# output to file
-#------------------------------------------------------------------------------
-
-start_screen_catcher(project.output_build_time_dir)
-
-print 'time stamp: ',project.gtime
-print 'USER: ', project.user
-
-
-copy_code_files(project.output_build_time_dir,__file__, 
-               dirname(project.__file__)+sep+ project.__name__+'.py' )
-
-
-
-#-------------------------------------------------------------------------------
-# 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
-#-------------------------------------------------------------------------------
-
-# topography directory filenames
-onshore_in_dir_name = project.onshore_in_dir_name
-coast_in_dir_name = project.coast_in_dir_name
-#island_in_dir_name = project.island_in_dir_name
-offshore_in_dir_name = project.offshore_in_dir_name
-
-onshore_dir_name = project.onshore_dir_name
-coast_dir_name = project.coast_dir_name
-#island_dir_name = project.island_dir_name
-offshore_dir_name = project.offshore_dir_name
-
-# creates DEM from asc data
-print "creates DEMs from ascii data"
-convert_dem_from_ascii2netcdf(onshore_in_dir_name, basename_out=onshore_dir_name, use_cache=True, verbose=True)
-#convert_dem_from_ascii2netcdf(island_in_dir_name, basename_out=island_dir_name, use_cache=True, verbose=True)
-
-#creates pts file for onshore DEM
-print "creates pts file for onshore DEM"
-dem2pts(onshore_dir_name,
-#        easting_min=project.eastingmin,
-#        easting_max=project.eastingmax,
-#        northing_min=project.northingmin,
-#        northing_max= project.northingmax,
-        use_cache=True, 
-        verbose=True)
-
-#creates pts file for island DEM
-#dem2pts(island_dir_name, use_cache=True, verbose=True)
-
-print'create Geospatial data1 objects from topographies'
-G1 = Geospatial_data(file_name = onshore_dir_name + '.pts')
-G2 = Geospatial_data(file_name = coast_in_dir_name + '.xya')
-#G3 = Geospatial_data(file_name = island_dir_name + '.pts')
-G_off = Geospatial_data(file_name = offshore_in_dir_name + '.xya')
-
-print'add all geospatial objects'
-G = G1 + G2 + G_off 
-
-print'clip combined geospatial object by bounding polygon'
-#G_clipped = G.clip(project.bounding_polygon)
-#FIXME: add a clip function to pts
-#print'shape of clipped data', G_clipped.get_data_points().shape
-
-print'export combined DEM file'
-if access(project.topographies_dir,F_OK) == 0:
-    mkdir (project.topographies_dir)
-G.export_points_file(project.combined_dir_name + '.xya')
-#G_clipped.export_points_file(project.combined_dir_name + '.xya')
-
-'''
-print'project.combined_dir_name + 1.xya',project.combined_dir_name + '1.xya'
-G_all=Geospatial_data(file_name = project.combined_dir_name + '1.xya')
-print'split'
-G_all_1, G_all_2 = G_all.split(.10)
-print'export 1'
-G_all_1.export_points_file(project.combined_dir_name+'_small1' + '.xya')
-print'export 2'
-G_all_2.export_points_file(project.combined_dir_name+'_other1' + '.xya')
-
-
-#-------------------------------------------------------------------------
-# Convert URS to SWW file for boundary conditions 
-#-------------------------------------------------------------------------
-print 'starting to create boundary conditions'
-boundaries_in_dir_name = project.boundaries_in_dir_name
-
-from anuga.shallow_water.data_manager import urs2sww
-
-print 'minlat=project.north_boundary, maxlat=project.south_boundary',project.north_boundary, project.south_boundary
-print 'minlon= project.west_boundary, maxlon=project.east_boundary',project.west_boundary, project.east_boundary
-
-#import sys; sys.exit()
-
-#if access(project.boundaries_dir,F_OK) == 0:
-#    mkdir (project.boundaries_dir)
-
-from caching import cache
-cache(urs2sww,
-      (boundaries_in_dir_name,
-       project.boundaries_dir_name1), 
-      {'verbose': True,
-       'minlat': project.south_boundary,
-       'maxlat': project.north_boundary,
-       'minlon': project.west_boundary,
-       'maxlon': project.east_boundary,
-#       'minlat': project.south,
-#       'maxlat': project.north,
-#       'minlon': project.west,
-#       'maxlon': project.east,
-       'mint': 0, 'maxt': 40000,
-#       'origin': domain.geo_reference.get_origin(),
-       'mean_stage': project.tide,
-#       'zscale': 1,                 #Enhance tsunami
-       'fail_on_NaN': False},
-       verbose = True,
-       )
-#       dependencies = source_dir + project.boundary_basename + '.sww')
-
-'''
-
-
-
-
-
-
-
Index: /anuga_work/production/geraldton/run_geraldton.py
===================================================================
--- /anuga_work/production/geraldton/run_geraldton.py	(revision 5007)
+++ /anuga_work/production/geraldton/run_geraldton.py	(revision 5007)
@@ -0,0 +1,296 @@
+"""Script for running tsunami inundation scenario for Dampier, 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.output_run_time_dir 
+
+The scenario is defined by a triangular mesh created from project.polygon,
+the elevation data and a simulated tsunami generated with URS code.
+
+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.shallow_water import Field_boundary
+from Numeric import allclose
+from anuga.shallow_water.data_manager import export_grid
+
+from anuga.pmesh.mesh_interface import create_mesh_from_regions
+from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
+from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
+from anuga_parallel.parallel_abstraction import get_processor_name
+from anuga.caching import myhash
+from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
+from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
+
+# Application specific imports
+import project                 # Definition of file names and polygons
+
+def run_model(**kwargs):
+    
+
+    #------------------------------------------------------------------------------
+    # Copy scripts to time stamped output directory and capture screen
+    # output to file
+    #------------------------------------------------------------------------------
+    print "Processor Name:",get_processor_name()
+
+    #copy script must be before screen_catcher
+    #print kwargs
+
+    print 'output_dir',kwargs['output_dir']
+    if myid == 0:
+        copy_code_files(kwargs['output_dir'],__file__, 
+                 dirname(project.__file__)+sep+ project.__name__+'.py' )
+
+        store_parameters(**kwargs)
+
+    barrier()
+
+    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
+
+    print "Processor Name:",get_processor_name()
+
+    # filenames
+#    meshes_dir_name = project.meshes_dir_name+'.msh'
+
+    # creates copy of code in output dir
+    print 'min triangles', project.trigs_min,
+    print 'Note: This is generally about 20% less than the final amount'
+
+    #--------------------------------------------------------------------------
+    # 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
+    #--------------------------------------------------------------------------
+
+    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
+    # causes problems with the ability to cache set quantity which takes alot of times
+    if myid == 0:
+    
+        print 'start create mesh from regions'
+
+        create_mesh_from_regions(project.poly_all,
+                             boundary_tags=project.boundary_tags,
+                             maximum_triangle_area=project.res_poly_all,
+                             interior_regions=project.interior_regions,
+                             filename=project.meshes_dir_name+'.msh',
+                             use_cache=False,
+                             verbose=True)
+    barrier()
+    
+
+    #-------------------------------------------------------------------------
+    # Setup computational domain
+    #-------------------------------------------------------------------------
+    print 'Setup computational domain'
+
+    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
+    #above don't work
+    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
+    print 'memory usage before del domain',mem_usage()
+       
+    print domain.statistics()
+    print 'triangles',len(domain)
+    
+    kwargs['act_num_trigs']=len(domain)
+
+
+    #-------------------------------------------------------------------------
+    # Setup initial conditions
+    #-------------------------------------------------------------------------
+    if myid == 0:
+
+        print 'Setup initial conditions'
+
+        from polygon import Polygon_function
+        #following sets the stage/water to be offcoast only
+        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
+                                 geo_reference = domain.geo_reference)
+        domain.set_quantity('stage', IC)
+#        domain.set_quantity('stage', kwargs['tide'])
+        domain.set_quantity('friction', kwargs['friction']) 
+        
+        print 'Start Set quantity',kwargs['bathy_file']
+
+        domain.set_quantity('elevation', 
+                            filename = kwargs['bathy_file'],
+                            use_cache = True,
+                            verbose = True,
+                            alpha = kwargs['alpha'])
+        print 'Finished Set quantity'
+    barrier()
+
+
+    #------------------------------------------------------
+    # Distribute domain to implement parallelism !!!
+    #------------------------------------------------------
+
+    if numprocs > 1:
+        domain=distribute(domain)
+
+    #------------------------------------------------------
+    # Set domain parameters
+    #------------------------------------------------------
+    print 'domain id', id(domain)
+    domain.set_name(kwargs['aa_scenario_name'])
+    domain.set_datadir(kwargs['output_dir'])
+    domain.set_default_order(2) # Apply second order scheme
+    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
+    domain.set_store_vertices_uniquely(False)
+    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
+    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    print 'domain id', id(domain)
+    domain.beta_h = 0
+    #domain.tight_slope_limiters = 1
+
+    #-------------------------------------------------------------------------
+    # Setup boundary conditions 
+    #-------------------------------------------------------------------------
+    print 'Available boundary tags', domain.get_boundary_tags()
+    print 'domain id', id(domain)
+    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
+
+    if project.source != 'other':
+        Bf = Field_boundary(kwargs['boundary_file'],
+                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
+                    use_cache=True, verbose=True)
+                    
+    kwargs['input_start_time']=domain.starttime
+
+    print 'finished reading boundary file'
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
+    Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
+
+
+    print'set_boundary'
+
+    domain.set_boundary({'back': Br,
+                         'side': Bd,
+                         'ocean': Bf}) 
+    print'finish set boundary'
+
+    #----------------------------------------------------------------------------
+    # Evolve system through time
+    #--------------------------------------------------------------------
+    t0 = time.time()
+
+    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
+        domain.write_time()
+        domain.write_boundary_statistics(tags = 'ocean')     
+
+        if allclose(t, 120):
+            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
+
+        if allclose(t, 720):
+            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
+
+#    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime']
+#                       ,skip_initial_step = True): 
+#        domain.write_time()
+#        domain.write_boundary_statistics(tags = 'ocean')     
+#    
+#    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
+#                       ,skip_initial_step = True): 
+#        domain.write_time()
+#        domain.write_boundary_statistics(tags = 'ocean')   
+
+    x, y = domain.get_maximum_inundation_location()
+    q = domain.get_maximum_inundation_elevation()
+
+    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
+
+    print 'That took %.2f seconds' %(time.time()-t0)
+
+    #kwargs 'completed' must be added to write the final parameters to file
+    kwargs['completed']=str(time.time()-t0)
+    
+    if myid==0:
+        store_parameters(**kwargs)
+    barrier
+    
+    print 'memory usage before del domain1',mem_usage()
+    
+def export_model(**kwargs):    
+    #store_parameters(**kwargs)
+    
+#    print 'memory usage before del domain',mem_usage()
+    #del domain
+    print 'memory usage after del domain',mem_usage()
+    
+    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
+    print'swwfile',swwfile
+    
+    export_grid(swwfile, extra_name_out = 'town',
+            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
+            timestep = None,
+            reduction = max,
+            cellsize = kwargs['export_cellsize'],
+            NODATA_value = -9999,
+            easting_min = project.eastingmin,
+            easting_max = project.eastingmax,
+            northing_min = project.northingmin,
+            northing_max = project.northingmax,
+            verbose = False,
+            origin = None,
+            datum = 'GDA94',
+            format = 'asc')
+            
+    inundation_damage(swwfile+'.sww', project.buildings_filename, 
+                      project.buildings_filename_out)
+    
+#-------------------------------------------------------------
+if __name__ == "__main__":
+    
+    kwargs={}
+    kwargs['est_num_trigs']=project.trigs_min
+    kwargs['num_cpu']=numprocs
+    kwargs['host']=project.host
+    kwargs['res_factor']=project.res_factor
+    kwargs['starttime']=project.starttime
+    kwargs['yieldstep']=project.yieldstep
+    kwargs['midtime']=project.midtime
+    kwargs['finaltime']=project.finaltime
+    
+    kwargs['output_dir']=project.output_run_time_dir
+    kwargs['bathy_file']=project.combined_dir_name
+#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
+    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
+    kwargs['file_name']=project.home+'detail.csv'
+    kwargs['aa_scenario_name']=project.scenario_name
+    kwargs['ab_time']=project.time
+    kwargs['res_factor']= project.res_factor
+    kwargs['tide']=project.tide
+    kwargs['user']=project.user
+    kwargs['alpha'] = project.alpha
+    kwargs['friction']=project.friction
+    kwargs['time_thinning'] = project.time_thinning
+    kwargs['dir_comment']=project.dir_comment
+    kwargs['export_cellsize']=project.export_cellsize
+    
+
+    run_model(**kwargs)
+      
+    if myid==0:
+        export_model(**kwargs)
+    barrier
Index: uga_work/production/geraldton/run_gerladton.py
===================================================================
--- /anuga_work/production/geraldton/run_gerladton.py	(revision 5006)
+++ 	(revision )
@@ -1,296 +1,0 @@
-"""Script for running tsunami inundation scenario for Dampier, 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.output_run_time_dir 
-
-The scenario is defined by a triangular mesh created from project.polygon,
-the elevation data and a simulated tsunami generated with URS code.
-
-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.shallow_water import Field_boundary
-from Numeric import allclose
-from anuga.shallow_water.data_manager import export_grid
-
-from anuga.pmesh.mesh_interface import create_mesh_from_regions
-from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
-from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
-from anuga_parallel.parallel_abstraction import get_processor_name
-from anuga.caching import myhash
-from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
-from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
-
-# Application specific imports
-import project                 # Definition of file names and polygons
-
-def run_model(**kwargs):
-    
-
-    #------------------------------------------------------------------------------
-    # Copy scripts to time stamped output directory and capture screen
-    # output to file
-    #------------------------------------------------------------------------------
-    print "Processor Name:",get_processor_name()
-
-    #copy script must be before screen_catcher
-    #print kwargs
-
-    print 'output_dir',kwargs['output_dir']
-    if myid == 0:
-        copy_code_files(kwargs['output_dir'],__file__, 
-                 dirname(project.__file__)+sep+ project.__name__+'.py' )
-
-        store_parameters(**kwargs)
-
-    barrier()
-
-    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
-
-    print "Processor Name:",get_processor_name()
-
-    # filenames
-#    meshes_dir_name = project.meshes_dir_name+'.msh'
-
-    # creates copy of code in output dir
-    print 'min triangles', project.trigs_min,
-    print 'Note: This is generally about 20% less than the final amount'
-
-    #--------------------------------------------------------------------------
-    # 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
-    #--------------------------------------------------------------------------
-
-    #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
-    # causes problems with the ability to cache set quantity which takes alot of times
-    if myid == 0:
-    
-        print 'start create mesh from regions'
-
-        create_mesh_from_regions(project.poly_all,
-                             boundary_tags=project.boundary_tags,
-                             maximum_triangle_area=project.res_poly_all,
-                             interior_regions=project.interior_regions,
-                             filename=project.meshes_dir_name+'.msh',
-                             use_cache=False,
-                             verbose=True)
-    barrier()
-    
-
-    #-------------------------------------------------------------------------
-    # Setup computational domain
-    #-------------------------------------------------------------------------
-    print 'Setup computational domain'
-
-    #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
-    #above don't work
-    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
-    print 'memory usage before del domain',mem_usage()
-       
-    print domain.statistics()
-    print 'triangles',len(domain)
-    
-    kwargs['act_num_trigs']=len(domain)
-
-
-    #-------------------------------------------------------------------------
-    # Setup initial conditions
-    #-------------------------------------------------------------------------
-    if myid == 0:
-
-        print 'Setup initial conditions'
-
-        from polygon import Polygon_function
-        #following sets the stage/water to be offcoast only
-        IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
-                                 geo_reference = domain.geo_reference)
-        domain.set_quantity('stage', IC)
-#        domain.set_quantity('stage', kwargs['tide'])
-        domain.set_quantity('friction', kwargs['friction']) 
-        
-        print 'Start Set quantity',kwargs['bathy_file']
-
-        domain.set_quantity('elevation', 
-                            filename = kwargs['bathy_file'],
-                            use_cache = True,
-                            verbose = True,
-                            alpha = kwargs['alpha'])
-        print 'Finished Set quantity'
-    barrier()
-
-
-    #------------------------------------------------------
-    # Distribute domain to implement parallelism !!!
-    #------------------------------------------------------
-
-    if numprocs > 1:
-        domain=distribute(domain)
-
-    #------------------------------------------------------
-    # Set domain parameters
-    #------------------------------------------------------
-    print 'domain id', id(domain)
-    domain.set_name(kwargs['aa_scenario_name'])
-    domain.set_datadir(kwargs['output_dir'])
-    domain.set_default_order(2) # Apply second order scheme
-    domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
-    domain.set_store_vertices_uniquely(False)
-    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
-    domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
-    print 'domain id', id(domain)
-    domain.beta_h = 0
-    #domain.tight_slope_limiters = 1
-
-    #-------------------------------------------------------------------------
-    # Setup boundary conditions 
-    #-------------------------------------------------------------------------
-    print 'Available boundary tags', domain.get_boundary_tags()
-    print 'domain id', id(domain)
-    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
-
-    if project.source != 'other':
-        Bf = Field_boundary(kwargs['boundary_file'],
-                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
-                    use_cache=True, verbose=True)
-                    
-    kwargs['input_start_time']=domain.starttime
-
-    print 'finished reading boundary file'
-
-    Br = Reflective_boundary(domain)
-    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
-    Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
-
-
-    print'set_boundary'
-
-    domain.set_boundary({'back': Br,
-                         'side': Bd,
-                         'ocean': Bf}) 
-    print'finish set boundary'
-
-    #----------------------------------------------------------------------------
-    # Evolve system through time
-    #--------------------------------------------------------------------
-    t0 = time.time()
-
-    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['starttime']): 
-        domain.write_time()
-        domain.write_boundary_statistics(tags = 'ocean')     
-
-        if allclose(t, 120):
-            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
-
-        if allclose(t, 720):
-            domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
-
-#    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime']
-#                       ,skip_initial_step = True): 
-#        domain.write_time()
-#        domain.write_boundary_statistics(tags = 'ocean')     
-#    
-#    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
-#                       ,skip_initial_step = True): 
-#        domain.write_time()
-#        domain.write_boundary_statistics(tags = 'ocean')   
-
-    x, y = domain.get_maximum_inundation_location()
-    q = domain.get_maximum_inundation_elevation()
-
-    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
-
-    print 'That took %.2f seconds' %(time.time()-t0)
-
-    #kwargs 'completed' must be added to write the final parameters to file
-    kwargs['completed']=str(time.time()-t0)
-    
-    if myid==0:
-        store_parameters(**kwargs)
-    barrier
-    
-    print 'memory usage before del domain1',mem_usage()
-    
-def export_model(**kwargs):    
-    #store_parameters(**kwargs)
-    
-#    print 'memory usage before del domain',mem_usage()
-    #del domain
-    print 'memory usage after del domain',mem_usage()
-    
-    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
-    print'swwfile',swwfile
-    
-    export_grid(swwfile, extra_name_out = 'town',
-            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
-            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
-            timestep = None,
-            reduction = max,
-            cellsize = kwargs['export_cellsize'],
-            NODATA_value = -9999,
-            easting_min = project.eastingmin,
-            easting_max = project.eastingmax,
-            northing_min = project.northingmin,
-            northing_max = project.northingmax,
-            verbose = False,
-            origin = None,
-            datum = 'GDA94',
-            format = 'asc')
-            
-    inundation_damage(swwfile+'.sww', project.buildings_filename, 
-                      project.buildings_filename_out)
-    
-#-------------------------------------------------------------
-if __name__ == "__main__":
-    
-    kwargs={}
-    kwargs['est_num_trigs']=project.trigs_min
-    kwargs['num_cpu']=numprocs
-    kwargs['host']=project.host
-    kwargs['res_factor']=project.res_factor
-    kwargs['starttime']=project.starttime
-    kwargs['yieldstep']=project.yieldstep
-    kwargs['midtime']=project.midtime
-    kwargs['finaltime']=project.finaltime
-    
-    kwargs['output_dir']=project.output_run_time_dir
-    kwargs['bathy_file']=project.combined_dir_name
-#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
-    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
-    kwargs['file_name']=project.home+'detail.csv'
-    kwargs['aa_scenario_name']=project.scenario_name
-    kwargs['ab_time']=project.time
-    kwargs['res_factor']= project.res_factor
-    kwargs['tide']=project.tide
-    kwargs['user']=project.user
-    kwargs['alpha'] = project.alpha
-    kwargs['friction']=project.friction
-    kwargs['time_thinning'] = project.time_thinning
-    kwargs['dir_comment']=project.dir_comment
-    kwargs['export_cellsize']=project.export_cellsize
-    
-
-    run_model(**kwargs)
-      
-    if myid==0:
-        export_model(**kwargs)
-    barrier
