Ignore:
Timestamp:
Sep 25, 2008, 2:24:38 PM (15 years ago)
Author:
kristy
Message:

Updated all scripts to reflect Perth format. Also added 250m scripts

File:
1 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/busselton/run_busselton.py

    r5669 r5786  
    1 """Script for running tsunami inundation scenario for Dampier, WA, Australia.
    2 
    3 Source data such as elevation and boundary data is assumed to be available in
    4 directories specified by project.py
    5 The output sww file is stored in project.output_run_time_dir
     1"""Script for running a tsunami inundation scenario for busselton, WA, Australia.
    62
    73The scenario is defined by a triangular mesh created from project.polygon,
    8 the elevation data and a simulated tsunami generated with URS code.
    9 
    10 Ole Nielsen and Duncan Gray, GA - 2005 and Jane Sexton, Nick Bartzis, GA - 2006
     4the elevation data is compiled into a pts file through build_busselton.py
     5and a simulated tsunami is generated through an sts file from build_boundary.py.
     6
     7Input: sts file (build_boundary.py for respective event)
     8       pts file (build_busselton.py)
     9       information from project file
     10Outputs: sww file stored in project.output_run_time_dir
     11The export_results_all.py and get_timeseries.py is reliant
     12on the outputs of this script
     13
     14Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
     15Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
    1116"""
    1217
    … …  
    1722# Standard modules
    1823from os import sep
     24import os
    1925from os.path import dirname, basename
    20 import os
    2126from os import mkdir, access, F_OK
    2227from shutil import copy
    … …  
    3237from Numeric import allclose
    3338from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
    34 
    3539from anuga.pmesh.mesh_interface import create_mesh_from_regions
    3640from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
    37 #from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
    3841from anuga_parallel.parallel_abstraction import get_processor_name
    3942from anuga.caching import myhash
    … …  
    4245from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
    4346from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
    44 
     47from polygon import Polygon_function
     48   
    4549# Application specific imports
    46 import project                 # Definition of file names and polygons
    47 
     50import project  # Definition of file names and polygons
    4851numprocs = 1
    4952myid = 0
    … …  
    5154def run_model(**kwargs):
    5255   
    53 
    5456    #------------------------------------------------------------------------------
    5557    # Copy scripts to time stamped output directory and capture screen
    … …  
    5961
    6062    #copy script must be before screen_catcher
    61     #print kwargs
    6263
    6364    print 'output_dir',kwargs['output_dir']
    64     if myid == 0:
    65         copy_code_files(kwargs['output_dir'],__file__,
    66                  dirname(project.__file__)+sep+ project.__name__+'.py' )
    67 
    68         store_parameters(**kwargs)
    69 
    70     #barrier()
     65   
     66    copy_code_files(kwargs['output_dir'],__file__,
     67             dirname(project.__file__)+sep+ project.__name__+'.py' )
     68
     69    store_parameters(**kwargs)
    7170
    7271    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
    7372
    7473    print "Processor Name:",get_processor_name()
    75 
     74   
    7675    #-----------------------------------------------------------------------
    7776    # Domain definitions
    … …  
    7978
    8079    # Read in boundary from ordered sts file
    81     urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
     80    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir_event,project.scenario_name))
    8281
    8382    # Reading the landward defined points, this incorporates the original clipping
    8483    # polygon minus the 100m contour
    85     landward_bounding_polygon = read_polygon(project.boundaries_dir+'landward_bounding_polygon.txt')
     84    landward_bounding_polygon = read_polygon(project.landward_dir)
    8685
    8786    # Combine sts polyline with landward points
    … …  
    9089    # counting segments
    9190    N = len(urs_bounding_polygon)-1
    92     boundary_tags={'back': [N+1,N+2,N+3,N+4, N+5], 'side': [N,N+6],'ocean': range(N)}
    93 
    94    
     91
     92    # boundary tags refer to project.landward 4 points equals 5 segments start at N
     93    boundary_tags={'back': [N+1,N+2,N+3,N+4, N+5], 'side': [N,N+6], 'ocean': range(N)}
     94
    9595    #--------------------------------------------------------------------------
    96     # Create the triangular mesh based on overall clipping polygon with a
    97     # tagged
     96    # Create the triangular mesh based on overall clipping polygon with a tagged
    9897    # boundary and interior regions defined in project.py along with
    9998    # resolutions (maximal area of per triangle) for each polygon
    10099    #--------------------------------------------------------------------------
    101100
    102     #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
     101    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
    103102    # causes problems with the ability to cache set quantity which takes alot of times
    104     if myid == 0:
    105    
    106         print 'start create mesh from regions'
    107103       
    108         create_mesh_from_regions(bounding_polygon,
    109                              boundary_tags=boundary_tags,
    110                              maximum_triangle_area=project.res_poly_all,
    111                              interior_regions=project.interior_regions,
    112                              filename=project.meshes_dir_name+'.msh',
    113                              use_cache=False,
    114                              verbose=False)
    115     #barrier()
    116    
    117 ##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file = kwargs['bathy_file'],
    118 ##                                    alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
    119 ##                                    mesh_file = project.meshes_dir_name+'.msh')
    120 ##        print 'optimal alpha', covariance_value,alpha       
    121 
     104    print 'start create mesh from regions'
     105
     106    create_mesh_from_regions(bounding_polygon,
     107                         boundary_tags=boundary_tags,
     108                         maximum_triangle_area=project.res_poly_all,
     109                         interior_regions=project.interior_regions,
     110                         filename=project.meshes_dir_name,
     111                         use_cache=True,
     112                         verbose=True)
     113   
    122114    #-------------------------------------------------------------------------
    123115    # Setup computational domain
    … …  
    125117    print 'Setup computational domain'
    126118
    127     domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
     119    domain = Domain(project.meshes_dir_name, use_cache=False, verbose=True)
    128120    print 'memory usage before del domain',mem_usage()
    129121       
    … …  
    137129    # Setup initial conditions
    138130    #-------------------------------------------------------------------------
    139     if myid == 0:
    140 
    141         print 'Setup initial conditions'
    142 
    143         from polygon import Polygon_function
    144         #following sets the stage/water to be offcoast only
    145         IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
    146                                  geo_reference = domain.geo_reference)
    147         domain.set_quantity('stage', IC)
    148 #        domain.set_quantity('stage', kwargs['tide'])
    149         domain.set_quantity('friction', kwargs['friction'])
    150        
    151         print 'Start Set quantity',kwargs['bathy_file']
    152 
    153         domain.set_quantity('elevation',
    154                             filename = kwargs['bathy_file'],
    155                             use_cache = True,
    156                             verbose = True,
    157                             alpha = kwargs['alpha'])
    158         print 'Finished Set quantity'
    159     #barrier()
    160 
    161 
    162     #------------------------------------------------------
    163     # Distribute domain to implement parallelism !!!
    164     #------------------------------------------------------
    165 
    166     if numprocs > 1:
    167         domain=distribute(domain)
     131    print 'Setup initial conditions'
     132
     133    # sets the initial stage in the offcoast region only
     134    IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
     135                             geo_reference = domain.geo_reference)
     136    domain.set_quantity('stage', IC)
     137    #domain.set_quantity('stage',kwargs['tide'] )
     138    domain.set_quantity('friction', kwargs['friction'])
     139   
     140    print 'Start Set quantity',kwargs['elevation_file']
     141
     142    domain.set_quantity('elevation',
     143                        filename = kwargs['elevation_file'],
     144                        use_cache = False,
     145                        verbose = True,
     146                        alpha = kwargs['alpha'])
     147    print 'Finished Set quantity'
     148
     149##   #------------------------------------------------------
     150##    # Distribute domain to implement parallelism !!!
     151##    #------------------------------------------------------
     152##
     153##    if numprocs > 1:
     154##        domain=distribute(domain)
    168155
    169156    #------------------------------------------------------
    … …  
    171158    #------------------------------------------------------
    172159    print 'domain id', id(domain)
    173     domain.set_name(kwargs['aa_scenario_name'])
     160    domain.set_name(kwargs['scenario_name'])
    174161    domain.set_datadir(kwargs['output_dir'])
    175     domain.set_default_order(2) # Apply second order scheme
    176     domain.set_minimum_storable_height(0.01) # Don't store anything less than 1cm
     162    domain.set_default_order(2)                 # Apply second order scheme
     163    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
    177164    domain.set_store_vertices_uniquely(False)
    178165    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    … …  
    186173    print 'domain id', id(domain)
    187174   
    188     boundary_urs_out=project.boundaries_dir_name
     175    boundary_urs_out=project.boundaries_dir_event + sep + project.scenario_name
    189176
    190177    Br = Reflective_boundary(domain)
    … …  
    192179   
    193180    print 'Available boundary tags', domain.get_boundary_tags()
    194     Bf = Field_boundary(boundary_urs_out+'.sts', # Change from file_boundary
    195                    domain, mean_stage=project.tide,
     181    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
     182                   domain, mean_stage= project.tide,
    196183                   time_thinning=1,
    197184                   default_boundary=Bd,
    … …  
    199186                   verbose = True,
    200187                   boundary_polygon=bounding_polygon)
    201    
    202     domain.set_boundary({'back': Bd,
     188
     189    domain.set_boundary({'back': Br,
    203190                         'side': Bd,
    204                          'ocean': Bf}) #changed from Bf to Bd for large wave
     191                         'ocean': Bf})
    205192
    206193    kwargs['input_start_time']=domain.starttime
    … …  
    214201
    215202    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
    216                             ,skip_initial_step = False ):
     203                       ,skip_initial_step = False):
    217204        domain.write_time()
    218205        domain.write_boundary_statistics(tags = 'ocean')
    219206
     207    # these outputs should be checked with the resultant inundation map
    220208    x, y = domain.get_maximum_inundation_location()
    221209    q = domain.get_maximum_inundation_elevation()
    222 
    223210    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)
    224211
    225     print 'That took %.2f seconds' %(time.time()-t0)
     212    print 'Simulation took %.2f seconds' %(time.time()-t0)
    226213
    227214    #kwargs 'completed' must be added to write the final parameters to file
    228215    kwargs['completed']=str(time.time()-t0)
    229    
    230     if myid==0:
    231         store_parameters(**kwargs)
    232     #barrier
    233    
     216     
     217    store_parameters(**kwargs)
     218     
    234219    print 'memory usage before del domain1',mem_usage()
    235 
    236  #-------------------------------------------------------------
     220   
     221   
     222#-------------------------------------------------------------
    237223if __name__ == "__main__":
    238224   
    239225    kwargs={}
    240     kwargs['est_num_trigs']=project.trigs_min
    241     kwargs['num_cpu']=numprocs
    242     kwargs['host']=project.host
    243     kwargs['res_factor']=project.res_factor
    244     kwargs['starttime']=project.starttime
    245     kwargs['yieldstep']=project.yieldstep
    246226    kwargs['finaltime']=project.finaltime
    247    
    248227    kwargs['output_dir']=project.output_run_time_dir
    249     kwargs['bathy_file']=project.combined_dir_name+'.txt'
    250     kwargs['file_name']=project.home+'detail.csv'
    251     kwargs['aa_scenario_name']=project.scenario_name
    252     kwargs['ab_time']=project.time
    253     kwargs['res_factor']= project.res_factor
     228    kwargs['elevation_file']=project.combined_dir_name+'.pts'
     229    kwargs['scenario_name']=project.scenario_name
    254230    kwargs['tide']=project.tide
    255     kwargs['user']=project.user
    256231    kwargs['alpha'] = project.alpha
    257232    kwargs['friction']=project.friction
    258     kwargs['time_thinning'] = project.time_thinning
    259     kwargs['dir_comment']=project.dir_comment
    260     kwargs['export_cellsize']=project.export_cellsize
    261    
    262 
     233    #kwargs['num_cpu']=numprocs
     234    #kwargs['host']=project.host
     235    #kwargs['starttime']=project.starttime
     236    #kwargs['yieldstep']=project.yieldstep
     237    #kwargs['user']=project.user
     238    #kwargs['time_thinning'] = project.time_thinning
     239     
    263240    run_model(**kwargs)
    264241     
    265     #barrier
    266    
     242   
Note: See TracChangeset for help on using the changeset viewer.