Changeset 2350


Ignore:
Timestamp:
Feb 7, 2006, 4:56:20 PM (18 years ago)
Author:
sexton
Message:

tidying up scripts for Sydney tsunami scenario

Location:
production/sydney_2006
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • production/sydney_2006/project.py

    r2317 r2350  
    4444#nminvix = 6283000
    4545
    46 basename = 'slump'
     46basename = 'slump2'
    4747
    4848if sys.platform == 'win32':
    49     home = '..\..\..\..'     #Sandpit's parent dir
     49    home = '..\..\..\..\..'     #Sandpit's parent dir
    5050else:   
    5151    home = expanduser('~')
     
    6969
    7070#Georeferencing
    71 from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees
     71#from pyvolution.coordinate_transforms.redfearn import degminsec2decimal_degrees
     72from coordinate_transforms.redfearn import degminsec2decimal_degrees
    7273
    7374#Origin of existing dem (FIXME: Temporary measure)
     
    7778x0_origin = 314036.58727982 #revised 100m and 25m data
    7879y0_origin = 6224951.2960092
    79 mesh_origin = (refzone, x0_origin, y0_origin)  # input from Neil's data
     80#mesh_origin = (refzone, x0_origin, y0_origin)  # input from Neil's data
    8081
    8182# define clipping polygon
  • production/sydney_2006/run_sydney_smf.py

    r2317 r2350  
    6868                                 project.combineddemname + '.pts')
    6969                                 
    70 # this allows the user to switch between different clipping polygons
    71 # print 'Which total zone are you interested in?'
    72 # mytest = int(raw_input('0 = all, 1 = harbour and 2 = botany bay     '))
    73 
    74 mytest = 0
    75 
    7670# Create Triangular Mesh
    7771# Overall clipping polygon and interior regions defined in project.py
    7872from pmesh.create_mesh import create_mesh_from_regions
    7973
    80 if mytest == 0:
    81     # for whole region
    82     interior_res = 5000
    83     interior_regions = [[project.harbour_polygon_2, interior_res],
    84                         [project.botanybay_polygon_2, interior_res]] # maximal area of per triangle
     74# for whole region
     75interior_res = 5000 # maximal area of per triangle
     76interior_regions = [[project.harbour_polygon_2, interior_res],
     77                    [project.botanybay_polygon_2, interior_res]]
    8578
    86     m = cache(create_mesh_from_regions,
    87               project.diffpolygonall,
    88               {'boundary_tags': {'bottom': [0],
    89                                  'right1': [1], 'right0': [2],
    90                                  'right2': [3], 'top': [4], 'left1': [5],
    91                                  'left2': [6], 'left3': [7]},
    92                'resolution': 100000,
    93                'filename': meshname,           
    94                'interior_regions': interior_regions},
    95               #evaluate=True,   
    96               verbose = True)
    97     #import sys; sys.exit()
    98    
    99     #Add elevation data to the mesh - this is in place of set_quantity
    100     mesh_elevname = 'elevation.msh'   
    101     cache(fit_to_mesh_file,(meshname,
    102                      #project.finedemname + '.pts',
    103                      #project.coarsedemname + '.pts',
    104                      project.combineddemname + '.pts',
    105                      mesh_elevname,
    106                      DEFAULT_ALPHA),
    107             {'verbose': True, 'expand_search': True, 'precrop': True},
    108           dependencies = [meshname],
    109           #evaluate = True,   
    110           verbose = False)
    111     meshname = mesh_elevname
    112    
     79m = cache(create_mesh_from_regions,
     80            project.diffpolygonall,
     81            {'boundary_tags': {'bottom': [0],
     82                            'right1': [1], 'right0': [2],
     83                            'right2': [3], 'top': [4], 'left1': [5],
     84                            'left2': [6], 'left3': [7]},
     85            'resolution': 100000,
     86            'filename': meshname,           
     87            'interior_regions': interior_regions},
     88            #evaluate=True,   
     89            verbose = True)
     90
     91#Add elevation data to the mesh - this is in place of set_quantity
     92#mesh_elevname = 'elevation.msh'   
     93#cache(fit_to_mesh_file,(meshname, project.combineddemname + '.pts',
     94#                        mesh_elevname, DEFAULT_ALPHA),
     95#      {'verbose': True, 'expand_search': True, 'precrop': True},
     96#      dependencies = [meshname],
     97#      #evaluate = True,   
     98#      verbose = False)
     99#meshname = mesh_elevname
     100
    113101# Setup domain
    114102domain = cache(pmesh_to_domain_instance, (meshname, Domain),
    115103               dependencies = [meshname],                     
    116                verbose = True)               
     104               verbose = True)
     105
     106# This section replaced with fit_to_mesh_file above
     107domain.set_quantity('elevation',
     108                     filename = project.combineddemname + '.pts',
     109#                     filename = project.finedemname + '.pts',
     110                     use_cache = True,
     111                     verbose = True)
     112
     113         
    117114       
    118115print 'Number of triangles = ', len(domain)
     
    122119domain.set_datadir(project.outputdir)
    123120domain.store = True
    124 #domain.quantities_to_be_stored = ['stage']
    125121domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    126122
    127 #print 'Do you want to run the slump scenario?'
    128 #q2 = int(raw_input('0 for yes, 1 for no     '))
    129 
    130 q2 = 0
    131 
    132 if q2 == 0:
    133     # slump parameters
    134     # as advised by Adrian, this can be as much as 15000-30000 metres
    135     # len = 4500.0
    136     len = 30000.0
    137     #thk = 670 #thk = 760.0
    138     dep = 400.0
    139     thk = 0.44*dep #thk = 120.0
    140     wid = 15000.0 # for slump scenario, wid=len
    141     #dep = 400 #1000.0 #d = 400 doesn't make sense with thk = 670!
    142     rad = 3330
    143     dp = 0.23
    144     th = 6
    145     alpha = 0.0
    146     x0 = project.x0
    147     y0 = project.y0
    148     # slide parameters
    149     # len = 10000.0
    150     # dep = 1200.0
    151     # thk = 200.0
    152     # th = 12
    153     # Slump or Slide Initial conditions
    154     t = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
    155                       radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
    156     #t = slide_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
    157     #                  x0=x0, y0=y0, alpha=alpha)
    158     domain.set_quantity('stage', t, verbose = True)
    159 
    160 if q2 == 1:
    161     domain.set_quantity('stage', tide)
    162 
    163 
    164 # These quantities are the same regardless of the initial condition
     123# slump parameters, wid=len
     124len = 30000.0
     125dep = 400.0
     126#thk = 120.0 for scenario 1 and 176 for scenario 2 rang 0.2-0.44 of depth
     127thk = 0.44*dep
     128rad = 3330
     129dp = 0.23
     130th = 6
     131alpha = 0.0
     132x0 = project.x0 # slump origin
     133y0 = project.y0
     134# Setup Initial conditions
     135t = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
     136                  radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
     137domain.set_quantity('stage', t)
    165138domain.set_quantity('friction', 0)
    166 
    167    
     139   
    168140# Setup Boundary Conditions
    169141print domain.get_boundary_tags()
     
    176148                   f=lambda t: [(6<t<606)*6, 0, 0])
    177149
    178 # for slump scenario
    179 if q2 == 0:
    180     domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
    181                                  'right2': Br, 'top': Br, 'left1': Br,
    182                                  'left2': Br, 'left3': Br} )
    183 
    184 # for square wave scenario
    185 if q2 == 1:
    186     domain.set_boundary( {'top': Br, 'bottom': Br, 'right': Bw, 'left': Br} )
     150domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
     151                      'right2': Br, 'top': Br, 'left1': Br,
     152                      'left2': Br, 'left3': Br} )
    187153
    188154# Evolve
     
    190156t0 = time.time()
    191157
    192 for t in domain.evolve(yieldstep = 120, finaltime = 18000): #120 10800 = 3hrs, 18000 = 5hrs
     158for t in domain.evolve(yieldstep = 120, finaltime = 18000):
    193159    domain.write_time()
    194     domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')    
     160    domain.write_boundary_statistics(tags = 'bottom')    
    195161   
    196162print 'That took %.2f seconds' %(time.time()-t0)
    197 
    198 #############################################################
    199 # testing system on a smaller mesh
    200     #m = create_mesh_from_regions(project.harbour_polygon,
    201     #                             boundary_tags={'bottom': [0], 'top': [2],
    202     #                                            'right': [1], 'left': [3]},
    203     #                             resolution=100000, filename=meshname)
    204    
    205 # This section replaced with fit_to_mesh_file above
    206 #domain.set_quantity('elevation',
    207 #                     #filename = project.combineddemname + '.pts',
    208 #                     filename = project.finedemname + '.pts',
    209                      #filename = project.datadir + 'bathy_land25_orig.pts',
    210 #                     use_cache = True,
    211 #                     verbose = True)
    212 
    213 
    214 # for switching between overall clipping regions - NOT USED
    215 
    216 if mytest == 1:
    217     # for harbour region
    218     south = project.hsouth
    219     north = project.hnorth
    220     west = project.hwest
    221     east = project.heast
    222 
    223     #interior_regions = [[project.harbour_polygon, 25000]] # maximal area of per triangle
    224     # meshname = project.meshname + 'harbour.msh'                   
    225     m = cache(create_mesh_from_regions,
    226               project.polygon_h,
    227               {'boundary_tags': {'bottom': [0], 'top': [2],
    228                                  'right': [1], 'left': [3]},
    229                'resolution': 50000,
    230                'filename': meshname},           
    231      #          'interior_regions': interior_regions},
    232               evaluate=True,
    233               verbose = True)
    234 
    235 if mytest == 2:
    236     # for botany bay region
    237     south = project.bsouth
    238     north = project.bnorth
    239     west = project.bwest
    240     east = project.beast
    241 
    242     #interior_regions = [[project.botanybay_polygon, 25000]] # maximal area of per triangle
    243     # meshname = project.meshname + 'bay.msh'                   
    244     m = cache(create_mesh_from_regions,
    245               project.polygon_bb,
    246               {'boundary_tags': {'bottom': [0], 'top': [2],
    247                                  'right': [1], 'left': [3]},
    248                'resolution': 50000,
    249                'filename': meshname},
    250      #          'interior_regions': interior_regions},
    251               evaluate=True,
    252               verbose = True)
    253    
Note: See TracChangeset for help on using the changeset viewer.