Ignore:
Timestamp:
Aug 13, 2008, 4:40:45 PM (17 years ago)
Author:
kristy
Message:

updated all Geraldton files to reflex the new boundary file (copied from Perth) and new bathymetry

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/production/geraldton/run_geraldton.py

    r5442 r5652  
    1717# Standard modules
    1818from os import sep
     19import os
    1920from os.path import dirname, basename
    2021from os import mkdir, access, F_OK
     
    3031from anuga.shallow_water import Field_boundary
    3132from Numeric import allclose
    32 from anuga.shallow_water.data_manager import export_grid
     33from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
    3334
    3435from anuga.pmesh.mesh_interface import create_mesh_from_regions
    3536from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
    36 from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
     37#from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
    3738from anuga_parallel.parallel_abstraction import get_processor_name
    3839from anuga.caching import myhash
    3940from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
    4041from anuga.fit_interpolate.benchmark_least_squares import mem_usage
     42from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
     43from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
     44from Scientific.IO.NetCDF import NetCDFFile
    4145
    4246# Application specific imports
    4347import project                 # Definition of file names and polygons
     48numprocs = 1
     49myid = 0
    4450
    4551def run_model(**kwargs):
     
    6268        store_parameters(**kwargs)
    6369
    64     barrier()
     70   # barrier()
    6571
    6672    start_screen_catcher(kwargs['output_dir'], myid, numprocs)
     
    6874    print "Processor Name:",get_processor_name()
    6975
    70     # filenames
    71 #    meshes_dir_name = project.meshes_dir_name+'.msh'
    72 
    73     # creates copy of code in output dir
    74     print 'min triangles', project.trigs_min,
    75     print 'Note: This is generally about 20% less than the final amount'
    76 
     76   
     77    #-----------------------------------------------------------------------
     78    # Domain definitions
     79    #-----------------------------------------------------------------------
     80
     81    # Read in boundary from ordered sts file
     82    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))
     83
     84    # Reading the landward defined points, this incorporates the original clipping
     85    # polygon minus the 100m contour
     86    landward_bounding_polygon = read_polygon(project.polygons_dir+'landward_boundary.txt')
     87
     88    # Combine sts polyline with landward points
     89    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
     90   
     91    # counting segments
     92    N = len(urs_bounding_polygon)-1
     93    boundary_tags={'back': [N+2,N+3], 'side': [N,N+1,N+4],'ocean': range(N)}
     94
     95   
    7796    #--------------------------------------------------------------------------
    7897    # Create the triangular mesh based on overall clipping polygon with a
     
    88107        print 'start create mesh from regions'
    89108
    90         create_mesh_from_regions(project.poly_all,
    91                              boundary_tags=project.boundary_tags,
     109        create_mesh_from_regions(bounding_polygon,
     110                             boundary_tags=boundary_tags,
    92111                             maximum_triangle_area=project.res_poly_all,
    93112                             interior_regions=project.interior_regions,
    94113                             filename=project.meshes_dir_name+'.msh',
    95                              use_cache=False,
     114                             use_cache=True,
    96115                             verbose=True)
    97     barrier()
    98    
     116   # barrier()
     117
     118##        covariance_value,alpha = find_optimal_smoothing_parameter (data_file= kwargs['elevation_file'],
     119##                                alpha_list=[0.001, 0.01, 0.1, 0.15, 0.2, 0.25, 0.3, 0.4, 0.5],
     120##                                mesh_file = project.meshes_dir_name+'.msh')
     121##        print 'optimal alpha', covariance_value,alpha       
    99122
    100123    #-------------------------------------------------------------------------
     
    103126    print 'Setup computational domain'
    104127
    105     #domain = cache(Domain, (meshes_dir_name), {'use_cache':True, 'verbose':True}, verbose=True)
    106     #above don't work
    107128    domain = Domain(project.meshes_dir_name+'.msh', use_cache=False, verbose=True)
    108129    print 'memory usage before del domain',mem_usage()
     
    123144        from polygon import Polygon_function
    124145        #following sets the stage/water to be offcoast only
    125         IC = Polygon_function( [(project.poly_mainland, -1.0)], default = kwargs['tide'],
     146        IC = Polygon_function( [(project.poly_mainland, 0)], default = kwargs['tide'],
    126147                                 geo_reference = domain.geo_reference)
    127148        domain.set_quantity('stage', IC)
    128 #        domain.set_quantity('stage', kwargs['tide'])
     149        #domain.set_quantity('stage',kwargs['tide'] )
    129150        domain.set_quantity('friction', kwargs['friction'])
    130151       
    131         print 'Start Set quantity',kwargs['bathy_file']
     152        print 'Start Set quantity',kwargs['elevation_file']
    132153
    133154        domain.set_quantity('elevation',
    134                             filename = kwargs['bathy_file'],
    135                             use_cache = True,
     155                            filename = kwargs['elevation_file'],
     156                            use_cache = False,
    136157                            verbose = True,
    137158                            alpha = kwargs['alpha'])
    138159        print 'Finished Set quantity'
    139     barrier()
    140 
     160    #barrier()
     161
     162##    #------------------------------------------------------
     163##    # Create x,y,z file of mesh vertex!!!
     164##    #------------------------------------------------------
     165##        coord = domain.get_vertex_coordinates()
     166##        depth = domain.get_quantity('elevation')
     167##       
     168##        # Write vertex coordinates to file
     169##        filename=project.vertex_filename
     170##        fid=open(filename,'w')
     171##        fid.write('x (m), y (m), z(m)\n')
     172##        for i in range(len(coord)):
     173##            pt=coord[i]
     174##            x=pt[0]
     175##            y=pt[1]
     176##            z=depth[i]
     177##            fid.write('%.6f,%.6f,%.6f\n' %(x, y, z))
     178##
    141179
    142180    #------------------------------------------------------
     
    157195    domain.set_store_vertices_uniquely(False)
    158196    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
     197    domain.tight_slope_limiters = 1
    159198    print 'domain id', id(domain)
    160 
    161199
    162200    #-------------------------------------------------------------------------
     
    165203    print 'Available boundary tags', domain.get_boundary_tags()
    166204    print 'domain id', id(domain)
    167     #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'
    168     print'set_boundary'
    169 
     205   
     206    boundary_urs_out=project.boundaries_dir_name
     207   
     208    print 'Available boundary tags', domain.get_boundary_tags()
     209    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
     210                   domain, mean_stage= project.tide,
     211                   time_thinning=1,
     212                   use_cache=True,
     213                   verbose = True,
     214                   boundary_polygon=bounding_polygon)
     215
     216   
    170217    Br = Reflective_boundary(domain)
    171218    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
    172     Bo = Dirichlet_boundary([kwargs['tide']+5.0,0,0])
    173 
    174     if project.source != 'other':
    175         Bf = Field_boundary(kwargs['boundary_file'],
    176                     domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
    177                     use_cache=True, verbose=True)
    178         print 'finished reading boundary file'
    179         domain.set_boundary({'back': Br,
    180                              'side': Bd,
    181                              'ocean': Bf})
    182     else:
    183         print 'set ocean'               
    184         domain.set_boundary({'back': Br,
    185                              'side': Bd,
    186                              'ocean': Bd})
    187 
    188 #    kwargs['input_start_time']=domain.starttime
    189 
    190 
     219
     220    fid = NetCDFFile(boundary_urs_out+'.sts', 'r')    #Open existing file for read
     221    sts_time=fid.variables['time'][:]+fid.starttime
     222    tmin=min(sts_time)
     223    tmax=max(sts_time)
     224    fid.close()
     225
     226    print 'Boundary end time ', tmax-tmin
     227   
     228##    Bf = Field_boundary(kwargs['boundary_file'],
     229##                domain, time_thinning=kwargs['time_thinning'], mean_stage=kwargs['tide'],
     230##                use_cache=False, verbose=True)
     231
     232    domain.set_boundary({'back': Br,
     233                         'side': Bd,
     234                         'ocean': Bf})
     235
     236    kwargs['input_start_time']=domain.starttime
    191237
    192238    print'finish set boundary'
     
    197243    t0 = time.time()
    198244
    199     for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']):
     245    for t in domain.evolve(yieldstep = project.yieldstep, finaltime = kwargs['finaltime']
     246                       ,skip_initial_step = False):
    200247        domain.write_time()
    201         domain.write_boundary_statistics(tags = 'ocean')     
    202 
    203         if allclose(t, 240):
    204             domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bo})
    205 
    206         if allclose(t, 1440):
    207             domain.set_boundary({'back': Br, 'side': Bd, 'ocean': Bd})
    208 
    209 #    for t in domain.evolve(yieldstep = kwargs['yieldstep'], finaltime = kwargs['midtime']
    210 #                       ,skip_initial_step = True):
    211 #        domain.write_time()
    212 #        domain.write_boundary_statistics(tags = 'ocean')     
    213 #   
    214 #    for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']
    215 #                       ,skip_initial_step = True):
    216 #        domain.write_time()
    217 #        domain.write_boundary_statistics(tags = 'ocean')   
    218 
     248        domain.write_boundary_statistics(tags = 'ocean')
     249
     250        if t >= tmax-tmin:
     251            print 'changed to tide boundary condition at ocean'
     252            domain.set_boundary({'ocean': Bd})
     253           
    219254    x, y = domain.get_maximum_inundation_location()
    220255    q = domain.get_maximum_inundation_elevation()
     
    229264    if myid==0:
    230265        store_parameters(**kwargs)
    231     barrier
     266   # barrier
    232267   
    233268    print 'memory usage before del domain1',mem_usage()
    234269   
    235 def export_model(**kwargs):   
    236     #store_parameters(**kwargs)
    237    
    238 #    print 'memory usage before del domain',mem_usage()
    239     #del domain
    240     print 'memory usage after del domain',mem_usage()
    241    
    242     swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
    243     print'swwfile',swwfile
    244    
    245     export_grid(swwfile, extra_name_out = 'town',
    246             quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    247             #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    248             timestep = None,
    249             reduction = max,
    250             cellsize = kwargs['export_cellsize'],
    251             NODATA_value = -9999,
    252             easting_min = project.eastingmin,
    253             easting_max = project.eastingmax,
    254             northing_min = project.northingmin,
    255             northing_max = project.northingmax,
    256             verbose = False,
    257             origin = None,
    258             datum = 'GDA94',
    259             format = 'asc')
    260            
    261     inundation_damage(swwfile+'.sww', project.buildings_filename,
    262                       project.buildings_filename_out)
    263270   
    264271#-------------------------------------------------------------
     
    272279    kwargs['starttime']=project.starttime
    273280    kwargs['yieldstep']=project.yieldstep
    274     kwargs['midtime']=project.midtime
    275281    kwargs['finaltime']=project.finaltime
    276282   
    277283    kwargs['output_dir']=project.output_run_time_dir
    278 #    kwargs['bathy_file']=project.combined_dir_name+'.txt'
    279     kwargs['bathy_file']=project.combined_dir_name_small + '.txt'
     284    kwargs['elevation_file']=project.combined_dir_name+'.pts'
    280285    kwargs['boundary_file']=project.boundaries_in_dir_name + '.sww'
    281286    kwargs['file_name']=project.home+'detail.csv'
     
    296301    if myid==0:
    297302        export_model(**kwargs)
    298     barrier
     303    #barrier
Note: See TracChangeset for help on using the changeset viewer.