Changeset 5111 for anuga_work/production

Timestamp:

Mar 4, 2008, 3:00:23 PM (17 years ago)

Author:

nick

Message:

updates to perth

Location:

anuga_work/production/perth

Files:

• project.py (modified) (6 diffs)
• run_perth.py (modified) (2 diffs)

anuga_work/production/perth/project.py

@@ -9 +9 @@
 from anuga.utilities.polygon import read_polygon, plot_polygons, is_inside_polygon, number_mesh_triangles
 #from anuga.coordinate_transforms.redfearn import degminsec2decimal_degrees, convert_points_from_latlon_to_utm
-from anuga.utilities.system_tools import get_user_name
+from anuga.utilities.system_tools import get_user_name, get_host_name
 
 # file and system info
 #---------------------------------
-codename = 'project.py'
-
-home = getenv('INUNDATIONHOME') #Sandpit's parent dir   
+#codename = 'project.py'
+
+home = getenv('INUNDATIONHOME') + sep +'data'+sep #Sandpit's parent diruser = get_user_name()
 user = get_user_name()
+host = get_host_name()
 
 # INUNDATIONHOME is the inundation directory, not the data directory.
-home += sep +'data'
 
 #time stuff
@@ -28 +28 @@
 print 'gtime: ', gtime
 
-tide = 0.6
-
 #Making assumptions about the location of scenario data
 state = 'western_australia'
 scenario_name = 'perth'
-scenario = 'perth_tsunami_scenario_2006'
+scenario = 'perth_tsunami_scenario'
+
+tide = 0.6
+
+alpha = 0.1
+friction=0.01
+starttime=10000
+midtime=21600
+finaltime=10000
+export_cellsize=50
+setup='trial'
+source='test'
+
+if setup =='trial':
+    print'trial'
+    res_factor=10
+    time_thinning=48
+    yieldstep=240
+if setup =='basic': 
+    print'basic'
+    res_factor=4
+    time_thinning=12
+    yieldstep=120
+if setup =='final': 
+    print'final'
+    res_factor=1
+    time_thinning=4
+    yieldstep=60
+
+dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user)
 
 # onshore data provided by WA DLI
@@ -52 +79 @@
 combined_smaller_name = 'perth_combined_elevation_smaller'
 
+anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
 
 topographies_in_dir = home+sep+state+sep+scenario+sep+'elevation_final'+sep+'points'+sep
-topographies_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'topographies'+sep
+topographies_dir = anuga_dir+'topographies'+sep
 topographies_time_dir = topographies_dir+build_time+sep
 
@@ -83 +111 @@
 #combined_time_dir_final_name = topographies_time_dir + combined_final_name
 
-
-
-meshes_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'meshes'+sep
+meshes_dir = anuga_dir+'meshes'+sep
 meshes_dir_name = meshes_dir + scenario_name
 
-polygons_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'polygons'+sep
-tide_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'tide_data'+sep
-
-
-boundaries_source = '1'
+polygons_dir = anuga_dir+'polygons'+sep
+tide_dir = anuga_dir+'tide_data'+sep
+
+
+#boundaries_source = '1'
+
+if source =='dampier':
+    boundaries_name = 'broome_3854_17042007' #Dampier gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'dampier'+sep+'1_10000'+sep
+
+if source=='onslow':
+    boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
+    
+if source=='exmouth':
+    boundaries_name = 'broome_3103_18052007' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep+'urs'+sep+'exmouth'+sep+'1_10000'+sep
+
+if source=='test':
+    boundaries_name = 'other' #exmouth gun
+    boundaries_in_dir = anuga_dir+'boundaries'+sep
+
+
 #boundaries locations
-boundaries_in_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+'urs'+sep+boundaries_source+sep
 boundaries_in_dir_name = boundaries_in_dir + scenario_name
-boundaries_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep
+boundaries_dir = anuga_dir+'boundaries'+sep
 boundaries_dir_name = boundaries_dir + scenario_name
-#boundaries_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'boundaries'+sep+build_time+sep
+#boundaries_time_dir = anuga_dir+'boundaries'+sep+build_time+sep
 #boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
 
 #output locations
-output_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep
-output_build_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+build_time+sep
-output_run_time_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'outputs'+sep+run_time+sep
+output_dir = anuga_dir+'outputs'+sep
+output_build_time_dir = anuga_dir+'outputs'+sep+build_time+sep
+output_run_time_dir = anuga_dir+'outputs'+sep+dir_comment+sep
 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
 
 #gauges
 gauge_name = 'perth.csv'
-gauges_dir = home+sep+state+sep+scenario+sep+'anuga'+sep+'gauges'+sep
+gauges_dir = anuga_dir+'gauges'+sep
 gauges_dir_name = gauges_dir + gauge_name
 
@@ -121 +164 @@
 
 poly_all = read_polygon(polygons_dir+'bounding_poly.csv')
-res_poly_all = 100000
-
-refzone = 50 
-
+res_poly_all = 100000*res_factor
+
+#refzone = 50 
 
 ###############################
@@ -132 +174 @@
 #poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20pts.csv')
 poly_pos20_neg20 = read_polygon(polygons_dir+'pos20_neg20_new_pts.csv')
-res_pos20_neg20 = 20000
+res_pos20_neg20 = 20000*res_factor
 
 #poly_cbd = read_polygon(polygons_dir+'cbd_pts.csv')
-poly_cbd = read_polygon(polygons_dir+'cbd_new_pts.csv')
-res_cbd = 1000
+poly_cbd = read_polygon(polygons_dir+'cbd_smaller_pts.csv')
+res_cbd = 500*res_factor
 
 poly_penguin = read_polygon(polygons_dir+'penguin_pts.csv')
-res_penguin = 1000
+res_penguin = 500*res_factor
+
+poly_geordie_bay = read_polygon(polygons_dir+'geordie_bay_pts.csv')
+res_geordie_bay = 500*res_factor
+
+poly_sorrento_gauge = read_polygon(polygons_dir+'sorrento_gauge_pts.csv')
+res_sorrento_gauge = 500*res_factor
 #assert zone == refzone
 
 interior_regions = [[poly_pos20_neg20,res_pos20_neg20],[poly_cbd,res_cbd]
+                     ,[poly_penguin,res_penguin],[poly_penguin,res_penguin]
                      ,[poly_penguin,res_penguin]]
+
+boundary_tags={'back': [4], 'side': [0,3],'ocean': [1, 2]}
 
 trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)

anuga_work/production/perth/run_perth.py

@@ -1 @@
-"""Script for running tsunami inundation scenario for Perth, WA, Australia.
+"""Script for running tsunami inundation scenario for Dampier, WA, Australia.
 
 Source data such as elevation and boundary data is assumed to be available in
@@ -28 +28 @@
 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.geospatial_data.geospatial_data import *
-from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
+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.utilities.polygon import plot_polygons, polygon_area
+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
 
-#------------------------------------------------------------------------------
-# Copy scripts to time stamped output directory and capture screen
-# output to file
-#------------------------------------------------------------------------------
-
-start_screen_catcher(project.output_run_time_dir, myid, numprocs)
-
-# filenames
-boundaries_name = project.scenario
-meshes_dir_name = project.meshes_dir_name+'.msh'
-#boundaries_time_dir_name = project.boundaries_time_dir_name
-boundaries_dir_name = project.boundaries_dir_name
-
-tide = project.tide
-print 'hello'
-# creates copy of code in output dir
-if myid == 0:
-    copy_code_files(project.output_run_time_dir,__file__, 
+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' )
-barrier()
-
-print 'USER: ', project.user
-print 'min triangles', project.trigs_min,
-print 'Note: This is generally about 20% less than the final amount'
-import sys
-sys.exit()
-#--------------------------------------------------------------------------
-# 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
-#--------------------------------------------------------------------------
-
-if myid == 0:
-
-    print 'start create mesh from regions'
-    create_mesh_from_regions(project.poly_all,
-                         boundary_tags={'back': [4], 'side': [0,3],
-                                        'ocean': [1, 2]},
-                         maximum_triangle_area=project.res_poly_all,
-                         interior_regions=project.interior_regions,
-#                         filename=meshes_time_dir_name,
-                         filename=meshes_dir_name,
-                         use_cache=True,
-                         verbose=True)
-
-# to sync all processors are ready 
-barrier()
-
-#-------------------------------------------------------------------------
-# Setup computational domain
-#-------------------------------------------------------------------------
-print 'Setup computational domain'
-#domain = Domain(meshes_time_dir_name, use_cache=True, verbose=True)
-domain = Domain(meshes_dir_name, use_cache=True, verbose=True)
-print domain.statistics()
-
-#-------------------------------------------------------------------------
-# Setup initial conditions
-#-------------------------------------------------------------------------
-if myid == 0:
-
-    print 'Setup initial conditions'
-
-    domain.set_quantity('stage', tide)
-    domain.set_quantity('friction', 0.01) 
-    #combined_time_dir_name = project.topographies_dir+build_time+project.combined_name
-    print 'Start Set quantity'
-
-    domain.set_quantity('elevation', 
-                    filename = project.combined_dir_name+'.txt',
-#                    filename = project.combined_smaller_name_dir+'.txt',
-                    use_cache = True,
-                    verbose = True,
-                    alpha = 0.1)
-    print 'Finished Set quantity'
-barrier()
-
-#------------------------------------------------------
-# Distribute domain to implement parallelism !!!
-#------------------------------------------------------
-
-if numprocs > 1:
-    domain=distribute(domain)
-
-#------------------------------------------------------
-# Set domain parameters
-#------------------------------------------------------
-
-domain.set_name(project.scenario_name)
-domain.set_datadir(project.output_run_time_dir)
-domain.set_default_order(2) # Apply second order scheme
-domain.set_minimum_storable_height(0.001) # Don't store anything less than 0.1cm
-#domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
-#domain.set_minimum_storable_height(0.1) # Don't store anything less than 10cm
-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)
-
-#-------------------------------------------------------------------------
-# Setup boundary conditions 
-#-------------------------------------------------------------------------
-print 'Available boundary tags', domain.get_boundary_tags()
-
-print 'Reading Boundary file'
-print 'domain id', id(domain)
-#Bf = File_boundary(boundaries_dir_name + '.sww',
-#                  domain, time_thinning=5, use_cache=True, verbose=True)
-
-print 'finished reading boundary file'
-Br = Reflective_boundary(domain)
-Bd = Dirichlet_boundary([tide,0,0])
-Bo = Dirichlet_boundary([tide+5.0,0,0])
-
-print'set_boundary'
-domain.set_boundary({'back': Br,
-                     'side': Bd,
-                     'ocean': Bd}) 
-print'finish set boundary'
-
-#----------------------------------------------------------------------------
-# Evolve system through time
-#----------------------------------------------------------------------------
-
-t0 = time.time()
-
-for t in domain.evolve(yieldstep = 60, finaltime = 10000):
-    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, 1020):
-        domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
-
-    
-print 'That took %.2f seconds' %(time.time()-t0)
-
+
+        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('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.tight_slope_limiters = 1
+    #domain.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
+    print 'domain id', id(domain)
+    domain.beta_h = 0
+
+    #-------------------------------------------------------------------------
+    # 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'
+
+    Br = Reflective_boundary(domain)
+    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
+    Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
+
+    if project.source != 'test':
+        print 'start reading boundary file'
+
+        Bf = Field_boundary(kwargs['boundary_file'],
+                    domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'], 
+                    use_cache=True, verbose=True)
+        domain.set_boundary({'back': Br,
+                             'side': Bd,
+                             'ocean': Bf}) 
+    else:
+        domain.set_boundary({'back': Br,
+                             'side': Bd,
+                             'ocean': Bd}) 
+
+    kwargs['input_start_time']=domain.starttime
+
+    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+'.txt'
+#    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