Changeset 5194


Ignore:
Timestamp:
Apr 9, 2008, 9:09:18 AM (16 years ago)
Author:
herve
Message:

run scripts for tonga

Location:
anuga_work/production/tonga
Files:
7 added
2 edited

Legend:

Unmodified
Added
Removed

  • anuga_work/production/tonga/project.py

    r5139 r5194  
    3030#Making assumptions about the location of scenario data
    3131state = 'sw_pacific'
    32 scenario_name = 'Tonga'
    33 scenario = 'fixed_wave'
     32scenario_name = 'Tonga_slide'
     33scenario = 'tonga'
    3434
    3535tide = 0
    … …  
    4242export_cellsize=50
    4343setup='trial'
     44source='test'
    4445
    4546if setup =='trial':
    … …  
    5960    yieldstep=60
    6061
    61 dir_comment='_'+setup+'_'+str(tide)+'_'+str(source)+'_'+str(user)
    62 
    63 
    64 topo_gridfile ='tongatapu_10mgrid'
     62dir_comment='_'+setup+'_'+str(tide)+'_'+str(scenario_name)+'_'+str(user)
     63
     64# onshore data 5m countour
     65onshore_name = 'topography_tongatapu' # original'
     66
     67# AHO + DPI data + colin French coastline
     68#coast_name = 'waterline'
     69Singlebeam_name = 'Tongatapu_SB_5m grid'
     70Multibeam_name = 'Tongatapu_MB_30m grid'
     71Chart_name= 'Tongatapu_Chart'
     72Derived_bath_name= 'Derived_Bathy'
     73added_data_name='joining_eastIsland_toreef'
     74
     75
     76#final topo name
     77combined_name ='Tongatapu_combined_elevation'
     78combined_smaller_name = 'Tongatapu_combined_elevation_smaller'
    6579
    6680anuga_dir = home+state+sep+scenario+sep+'anuga'+sep
    … …  
    7084#topographies_time_dir = topographies_dir+build_time+sep
    7185
     86# input topo file location
     87onshore_in_dir_name = topographies_in_dir + onshore_name
     88
     89Singlebeam_in_dir_name = topographies_in_dir + Singlebeam_name
     90Multibeam_in_dir_name = topographies_in_dir + Multibeam_name
     91Chart_in_dir_name = topographies_in_dir + Chart_name
     92Derived_bath_in_dir_name = topographies_in_dir + Derived_bath_name
     93added_data_in_dir_name = topographies_in_dir + added_data_name
     94
     95onshore_dir_name = topographies_dir + onshore_name
     96
     97Singlebeam_dir_name = topographies_dir + Singlebeam_name
     98Multibeam_dir_name = topographies_dir + Multibeam_name
     99Chart_dir_name = topographies_dir + Chart_name
     100Derived_bath_dir_name = topographies_dir + Derived_bath_name
     101added_data_dir_name = topographies_dir + added_data_name
     102
     103#final topo files
     104combined_dir_name = topographies_dir + combined_name
     105#combined_time_dir_name = topographies_time_dir + combined_name
     106combined_smaller_name_dir = topographies_dir + combined_smaller_name
     107#combined_time_dir_final_name = topographies_time_dir + combined_final_name
     108
    72109meshes_dir = anuga_dir+'meshes'+sep
    73110meshes_dir_name = meshes_dir + scenario_name
    … …  
    76113tide_dir = anuga_dir+'tide_data'+sep
    77114
     115
     116#boundaries_source = '1'
     117
     118if source =='dampier':
     119    boundaries_name = 'broome_3854_17042007' #Dampier gun
     120    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'dampier'+sep+'1_10000'+sep
     121
     122if source=='onslow':
     123    boundaries_name = 'broome_3859_16052007' #onslow_hedland_broome gun
     124    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'onslow_hedland_broome'+sep+'1_10000'+sep
     125   
     126if source=='exmouth':
     127    boundaries_name = 'broome_3103_18052007' #exmouth gun
     128    boundaries_in_dir = anuga_dir+'boundaries'+sep+sep+'exmouth'+sep+'1_10000'+sep
     129
     130if source=='test':
     131    boundaries_name = 'other' #exmouth gun
     132    boundaries_in_dir = anuga_dir+'boundaries'+sep
     133
     134
     135#boundaries locations
    78136boundaries_in_dir_name = boundaries_in_dir + scenario_name
    79137boundaries_dir = anuga_dir+'boundaries'+sep
    80138boundaries_dir_name = boundaries_dir + scenario_name
     139#boundaries_time_dir = anuga_dir+'boundaries'+sep+build_time+sep
     140#boundaries_time_dir_name = boundaries_time_dir + boundaries_name  #Used by post processing
    81141
    82142#output locations
    … …  
    84144output_build_time_dir = anuga_dir+'outputs'+sep+build_time+dir_comment+sep
    85145output_run_time_dir = anuga_dir+'outputs'+sep+run_time+dir_comment+sep
    86 output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post pr
     146output_run_time_dir_name = output_run_time_dir + scenario_name  #Used by post processing
     147
     148#gauges
     149#gauge_name = 'perth.csv'
     150#gauges_dir = anuga_dir+'gauges'+sep
     151#gauges_dir_name = gauges_dir + gauge_name
     152
     153#buildings_filename = gauges_dir + 'Perth_res_Project.csv'
     154#buildings_filename_out = 'Perth_res_Project_modified.csv'
    87155
    88156###############################
    … …  
    92160
    93161poly_all = read_polygon(polygons_dir+'extent.txt')
    94 res_poly_all = 100000*res_factor
    95 
     162res_poly_all = 40000*res_factor
     163
     164#refzone = 50
    96165
    97166###############################
    … …  
    99168###############################
    100169
    101 # interior polygons
    102 poly_tongatapu = read_polygon('tongatapu_finemesharea.txt')
    103 poly_island1 = read_polygon('island1.txt')
    104 poly_island2 = read_polygon('island2.txt')
    105 poly_island3 = read_polygon('islands3.txt')
    106 poly_island4 = read_polygon('islands4.txt')
    107 poly_island5= read_polygon('islands5.txt')
    108 poly_island6= read_polygon('islands6.txt')
    109 poly_island7= read_polygon('islands7.txt')
    110 poly_island8= read_polygon('islands8.txt')
    111 poly_island9= read_polygon('islands9.txt')
    112 poly_island10= read_polygon('islands10.txt')
    113 poly_island11= read_polygon('islands11.txt')
    114 poly_island12= read_polygon('islands12.txt')
    115 poly_island13= read_polygon('islands13.txt')
     170poly_tongatapu = read_polygon(polygons_dir+'finres.txt')
     171poly_island1 = read_polygon(polygons_dir+'island1.txt')
     172poly_island2 = read_polygon(polygons_dir+'island2.txt')
     173poly_island3 = read_polygon(polygons_dir+'island3.txt')
     174poly_island4 = read_polygon(polygons_dir+'island4.txt')
     175poly_island5= read_polygon(polygons_dir+'island5.txt')
     176poly_island6= read_polygon(polygons_dir+'island6.txt')
     177poly_island7= read_polygon(polygons_dir+'island7.txt')
     178poly_island8= read_polygon(polygons_dir+'island8.txt')
     179poly_island9= read_polygon(polygons_dir+'island9.txt')
     180poly_island10= read_polygon(polygons_dir+'island10.txt')
     181poly_island11= read_polygon(polygons_dir+'island11.txt')
     182poly_island12= read_polygon(polygons_dir+'island12.txt')
     183poly_island13= read_polygon(polygons_dir+'island13.txt')
    116184
    117185
    118186islands_res = 5000
    119 cairns_res = 5000
    120 interior_regions = [[project.poly_tongatapu, tongatapu_res],
    121                     [project.poly_island1, islands_res],
    122                     [project.poly_island2, islands_res],
    123                     [project.poly_island3, islands_res],
    124                     [project.poly_island4, islands_res],
    125                     [project.poly_island5, islands_res],
    126                     [project.poly_island6, islands_res],
    127                     [project.poly_island7, islands_res],
    128                     [project.poly_island8, islands_res],
    129                     [project.poly_island9, islands_res],
    130                     [project.poly_island10, islands_res],
    131                     [project.poly_island11, islands_res],
    132                     [project.poly_island12, islands_res],
    133                     [project.poly_island13, islands_res]]
    134 
    135 boundary_tags={'ocean_east': [0],'ocean_north': [1],'ocean_west': [2],'land': [3,4,5,6,7],'ocean_southeast':[8,9]}
     187tongatapu_res = 5000
     188interior_regions = [[poly_tongatapu, tongatapu_res],
     189                    [poly_island1, islands_res],
     190                    [poly_island2, islands_res],
     191                    [poly_island3, islands_res],
     192                    [poly_island4, islands_res],
     193                    [poly_island5, islands_res],
     194                    [poly_island6, islands_res],
     195                    [poly_island7, islands_res],
     196                    [poly_island8, islands_res],
     197                    [poly_island9, islands_res],
     198                    [poly_island10, islands_res],
     199                    [poly_island11, islands_res],
     200                    [poly_island12, islands_res],
     201                    [poly_island13, islands_res]]
     202
     203boundary_tags={'ocean_east':[0],
     204               'ocean_north':[1],
     205               'ocean_west':[2,3],
     206               'land':[4,5,6,7,8,9,10,11,12,13,14],
     207               'ocean_southeast':[15,16]}
    136208
    137209trigs_min = number_mesh_triangles(interior_regions, poly_all, res_poly_all)
    138210
    139211print 'min number triangles', trigs_min
    140 
    141212###################################################################
    142213# Clipping regions for export to asc and regions for clipping data
    … …  
    156227
    157228
    158 
    159 
    160 
    161 
    162 
    163 
    164 
    165 # topo Filenames
    166 dem_name = 'tongatapu_10mgrid'
    167 meshname = 'tongatapu.msh'
    168 
    169 # Create DEM from asc data
    170 convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
    171 
    172 # Create pts file for onshore DEM
    173 dem2pts(dem_name, use_cache=True, verbose=True)
    174 
    175 
    176 
    177 
    178 
    179 
    180 # -*- coding: cp1252 -*-
    181 """Common filenames and locations for topographic data, meshes and outputs.
    182 """
    183 
    184 from os import sep, environ, getenv, getcwd
    185 from os.path import expanduser
    186 import sys
    187 from time import localtime, strftime, gmtime
    188 from anuga.utilities.polygon import read_polygon, plot_polygons, \
    189                                     polygon_area, is_inside_polygon
    190 
    191 ###############################
    192 # Domain definitions
    193 ###############################
    194 
    195 # bounding polygon for study area
    196 bounding_polygon = read_polygon('Y:\data\sw_pacific\tonga\anuga\boundaries\extent.txt')
    197 
    198 print 'Area of bounding polygon', polygon_area(bounding_polygon)/1000000.0
    199 
    200 ###############################
    201 # Interior region definitions
    202 ###############################
    203 
    204 # interior polygons
    205 poly_tongatapu = read_polygon('tongatapu_finemesharea.txt')
    206 poly_island1 = read_polygon('island1.txt')
    207 poly_island2 = read_polygon('island2.txt')
    208 poly_island3 = read_polygon('islands3.txt')
    209 poly_island4 = read_polygon('islands4.txt')
    210 poly_island5= read_polygon('islands5.txt')
    211 poly_island6= read_polygon('islands6.txt')
    212 poly_island7= read_polygon('islands7.txt')
    213 poly_island8= read_polygon('islands8.txt')
    214 poly_island9= read_polygon('islands9.txt')
    215 poly_island10= read_polygon('islands10.txt')
    216 poly_island11= read_polygon('islands11.txt')
    217 poly_island12= read_polygon('islands12.txt')
    218 poly_island13= read_polygon('islands13.txt')
    219 
    220 #plot_polygons([bounding_polygon,poly_cairns,poly_island0,poly_island1,\
    221 #               poly_island2,poly_island3,poly_shallow],\
    222 #              'boundingpoly',verbose=False)
    223 
    224 ###################################################################
    225 # Clipping regions for export to asc and regions for clipping data
    226 ###################################################################
    227 
    228 # exporting asc grid
    229 eastingmin = 670500
    230 eastingmax = 712750
    231 northingmin = 7646000
    232 northingmax = 7677000
    233 
    234 
    235 slide_origin = [701290, 7665750] # move onto the continental shelf, depth = 500
    236 slide_depth = 207.
    237 
    238 #gauge_filename = 'gauges.csv'

  • anuga_work/production/tonga/tongatapu.py

    r5139 r5194  
    2424
    2525# Related major packages
     26from anuga.shallow_water import Time_boundary
    2627from anuga.shallow_water import Domain
    2728from anuga.shallow_water import Dirichlet_boundary
    … …  
    3132from Numeric import allclose
    3233from anuga.shallow_water.data_manager import export_grid
    33 
     34from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
    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
    … …  
    8182    # resolutions (maximal area of per triangle) for each polygon
    8283    #--------------------------------------------------------------------------
    83 
    84     #IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
    85     # causes problems with the ability to cache set quantity which takes alot of times
    8684   
    8785    if myid == 0:
    … …  
    9896    barrier()
    9997
     98    scenario='fixed_wave'
     99   
    100100    #-------------------------------------------------------------------------
    101101    # Setup computational domain
    … …  
    133133        domain.set_quantity('elevation',
    134134                            filename = kwargs['bathy_file'],
    135                             use_cache = True,
     135                            use_cache = False,
    136136                            verbose = True,
    137137                            alpha = kwargs['alpha'])
    … …  
    170170        tsunami_source = slide_tsunami(length=35000.0,
    171171                                       depth=project.slide_depth,
    172                                        slope=6.0,
    173                                        thickness=500.0,
     172                                       slope=10.0,
     173                                       thickness=800.0,
    174174                                       x0=project.slide_origin[0],
    175175                                       y0=project.slide_origin[1],
    … …  
    192192    # FIXME (Ole): Change this to Transmissive Momentum as
    193193    # suggested by Rajaraman (ticket:208)
     194        from math import sin, pi
    194195        Bw = Time_boundary(domain = domain,
    195                            f=lambda t: [(60<t<3660)*100, 0, 0])
    196         domain.set_boundary({'ocean_east': Bw,
    197                              'ocean_southeast': Bd,
    198                              'land': Bd,
    199                              'ocean_north': Bd,
    200                              'ocean_west':})
     196                           f=lambda t: [5.0*sin(t*pi*1/15), 0, 0])
     197        domain.set_boundary({'ocean_east':Bw,
     198                             'ocean_north':Bd,
     199                             'ocean_west':Bd,
     200                             'land':Bd,
     201                             'ocean_southeast':Bd})
    201202
    202203    # boundary conditions for slide scenario
    203204    if scenario == 'slide':
    204         domain.set_boundary({'ocean_east': Bd,
    205                              'ocean_southeast': Bd,
    206                              'land': Bd,
    207                              'ocean_north': Bd,
    208                              'ocean_west':})
     205        domain.set_boundary({'ocean_east':Bd,
     206                             'ocean_north':Bd,
     207                             'ocean_west':Bd,
     208                             'land':Bd,
     209                             'ocean_southeast':Bd})
    209210
    210211    #----------------------------------------------------------------------------
    … …  
    213214    t0 = time.time()
    214215
    215     for t in domain.evolve(yieldstep = 240, finaltime = kwargs['finaltime']):
     216    for t in domain.evolve(yieldstep = 1, finaltime = kwargs['finaltime']):
    216217        domain.write_time()
    217218        domain.write_boundary_statistics(tags = 'ocean_east')     
    … …  
    240241    print 'memory usage after del domain',mem_usage()
    241242   
    242     swwfile = kwargs['output_dir']+kwargs['scenario_name']
     243    swwfile = kwargs['output_dir']+kwargs['aa_scenario_name']
    243244    print'swwfile',swwfile
    244245   
    245     export_grid(swwfile, extra_name_out = 'town',
     246    export_grid(swwfile, extra_name_out = 'tonga',
    246247            quantities = ['speed','depth','elevation','stage'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    247248            #quantities = ['speed','depth'], # '(xmomentum**2 + ymomentum**2)**0.5' defaults to elevation
    248249            timestep = None,
    249250            reduction = max,
    250             cellsize = kwargs['export_cellsize'],
     251            cellsize = 25,
    251252            NODATA_value = -1E-030,
    252253            easting_min = project.eastingmin,
    … …  
    271272    kwargs['midtime']=project.midtime
    272273    kwargs['finaltime']=project.finaltime
    273    
    274274    kwargs['output_dir']=project.output_run_time_dir
    275275    kwargs['bathy_file']=project.combined_dir_name+'.txt'
Note: See TracChangeset for help on using the changeset viewer.