Changeset 2292


Ignore:
Timestamp:
Jan 26, 2006, 10:40:35 AM (18 years ago)
Author:
sexton
Message:

version 25 Jan - slump initial condition, new data

File:
1 edited

Legend:

Unmodified
Added
Removed
  • production/sydney_2006/run_sydney_smf.py

    r2240 r2292  
    2121from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
    2222     dem2pts
    23 from pyvolution.pmesh2domain import pmesh_to_domain_instance
     23from pyvolution.pmesh2domain import pmesh_to_domain_instance
     24from pyvolution.combine_pts import combine_rectangular_points_files
    2425from caching import cache
    2526import project
    2627from math import pi, sin
    2728from smf import slump_tsunami, slide_tsunami, Double_gaussian
    28 
    29 #Data preparation
    30 #Convert ASC 2 DEM 2 PTS using source data and store result in source data
    31 demname = project.demname
     29from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA
     30
     31# Data preparation
     32# Convert ASC 2 DEM 2 PTS using source data and store result in source data
     33# Do for coarse and fine data
     34# Fine pts file to be clipped to area of interest
     35coarsedemname = project.coarsedemname
     36finedemname = project.finedemname
    3237meshname = project.meshname+'.msh'
    3338
    34 cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
    35       dependencies = [demname + '.asc'],
     39# coarse data
     40cache(convert_dem_from_ascii2netcdf, coarsedemname, {'verbose': True},
     41      dependencies = [coarsedemname + '.asc'],
    3642      verbose = True)
    3743      #evaluate = True)
    3844
    39 cache(dem2pts, demname, {'verbose': True},
    40       dependencies = [demname + '.dem'],     
    41       verbose = True)
    42 
    43 #this allows the user to switch between different clipping polygons
    44 #print 'Which total zone are you interested in?'
    45 #mytest = int(raw_input('0 = all, 1 = harbour and 2 = botany bay     '))
     45cache(dem2pts, coarsedemname, {'verbose': True},
     46      dependencies = [coarsedemname + '.dem'],     
     47      verbose = True)
     48
     49# fine data
     50cache(convert_dem_from_ascii2netcdf, finedemname, {'verbose': True},
     51      dependencies = [finedemname + '.asc'],
     52      verbose = True)
     53      #evaluate = True)
     54
     55# clipping the fine data
     56cache(dem2pts, finedemname, {'verbose': True,
     57      'easting_min': project.eastingmin,
     58      'easting_max': project.eastingmax,
     59      'northing_min': project.northingmin,
     60      'northing_max': project.northingmax},
     61      dependencies = [finedemname + '.dem'],
     62      #evaluate = True,
     63      verbose = True)
     64
     65# combining the coarse and fine data
     66combine_rectangular_points_files(project.finedemname + '.pts',
     67                                 project.coarsedemname + '.pts',
     68                                 project.combineddemname + '.pts')
     69                                 
     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     '))
    4673
    4774mytest = 0
    4875
    49 #Create Triangular Mesh
     76# Create Triangular Mesh
     77# Overall clipping polygon and interior regions defined in project.py
    5078from pmesh.create_mesh import create_mesh_from_regions
    5179
    5280if mytest == 0:
    5381    # for whole region
    54     #south = project.south
    55     #north = project.north
    56     #west = project.west
    57     #east = project.east
    58 
    59     interior_regions = [[project.harbour_polygon_2, 25000],
    60                         [project.botanybay_polygon_2, 25000]] # maximal area of per triangle
    61 
    62     # meshname = project.meshname + 'all.msh'
     82    interior_res = 2000 # run at 2000 for final one
     83    interior_regions = [[project.harbour_polygon_2, interior_res],
     84                        [project.botanybay_polygon_2, interior_res]] # maximal area of per triangle
    6385
    6486    m = cache(create_mesh_from_regions,
    65               #project.polygonall,
    6687              project.diffpolygonall,
    67               {'boundary_tags': {'bottom': [0], 'top': [4],
    68                                  'right': [1], 'left1': [5],
    69                                  'left2': [6], 'left3': [7],
     88              {'boundary_tags': {'bottom': [0],
    7089                                 'right1': [1], 'right0': [2],
    71                                  'right2': [3]},
     90                                 'right2': [3], 'top': [4], 'left1': [5],
     91                                 'left2': [6], 'left3': [7]},
    7292               'resolution': 100000,
    7393               'filename': meshname,           
     
    7696              verbose = True)
    7797    #import sys; sys.exit()
    78 
    79 #Setup domain
    80 
     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   
     113# Setup domain
    81114domain = cache(pmesh_to_domain_instance, (meshname, Domain),
    82115               dependencies = [meshname],                     
     
    86119print 'The extent is ', domain.get_extent()
    87120
    88 #outputfilename = basename + '_slump'
    89 
    90121domain.set_name(project.basename)
    91 #domain.set_name(project.outfilename)
    92122domain.set_datadir(project.outputdir)
    93123domain.store = True
    94 
    95 #domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    96 domain.quantities_to_be_stored = ['stage']
     124#domain.quantities_to_be_stored = ['stage']
     125domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
    97126
    98127#print 'Do you want to run the slump scenario?'
     
    102131
    103132if q2 == 0:
    104     #slump parameters
    105     #as advised by Adrian, this can be as much as 15000-30000 metres
     133    # slump parameters
     134    # as advised by Adrian, this can be as much as 15000-30000 metres
    106135    # len = 4500.0
    107136    len = 15000.0
    108137    thk = 670 #thk = 760.0
    109     #wid = 4500.0
    110138    wid = 15000.0 # for slump scenario, wid=len
    111     dep = 1000.0 #d = 1200
     139    #dep = 400 #1000.0 #d = 400 doesn't make sense with thk = 670!
     140    dep = 1000.0
    112141    rad = 3330
    113142    dp = 0.23
    114143    th = 6
    115 
     144    alpha = 0.0
     145    x0 = project.x0
     146    y0 = project.y0
    116147    # slide parameters
    117148    # len = 10000.0
     
    119150    # thk = 200.0
    120151    # th = 12
    121     x0 = project.x0
    122     y0 = project.y0
    123     alpha = 0.0
    124     ######################
    125     # Slump Initial conditions
     152    # Slump or Slide Initial conditions
    126153    t = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
    127154                      radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
    128155    #t = slide_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
    129156    #                  x0=x0, y0=y0, alpha=alpha)
    130     domain.set_quantity('stage', t)
     157    domain.set_quantity('stage', t, verbose = True)
    131158
    132159if q2 == 1:
     
    136163# These quantities are the same regardless of the initial condition
    137164domain.set_quantity('friction', 0)
    138 domain.set_quantity('elevation',
    139                      filename = demname + '.pts',
    140                      use_cache = True,
    141                      verbose = True)
     165
    142166   
    143167# Setup Boundary Conditions
     
    151175                   f=lambda t: [(6<t<606)*6, 0, 0])
    152176
     177# for slump scenario
     178if q2 == 0:
     179    domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
     180                                 'right2': Br, 'top': Br, 'left1': Br,
     181                                 'left2': Br, 'left3': Br} )
     182
    153183# for square wave scenario
    154184if q2 == 1:
    155185    domain.set_boundary( {'top': Br, 'bottom': Br, 'right': Bw, 'left': Br} )
    156186
    157 # for slump scenario
    158 if q2 == 0:
    159     domain.set_boundary( {'top': Br, 'bottom': Br, 'right0': Br,
    160                           'left1': Br, 'left2': Br, 'left3': Br,
    161                           'right1': Br, 'right2': Br} )
    162 
    163187# Evolve
    164188import time
    165189t0 = time.time()
    166190
    167 for t in domain.evolve(yieldstep = 100, finaltime = 10400):
     191for t in domain.evolve(yieldstep = 120, finaltime = 10800): #120 10800 = 3hrs
    168192    domain.write_time()
    169193    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')     
    170 #    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')
    171194   
    172195print 'That took %.2f seconds' %(time.time()-t0)
     196
     197#############################################################
     198# testing system on a smaller mesh
     199    #m = create_mesh_from_regions(project.harbour_polygon,
     200    #                             boundary_tags={'bottom': [0], 'top': [2],
     201    #                                            'right': [1], 'left': [3]},
     202    #                             resolution=100000, filename=meshname)
     203   
     204# This section replaced with fit_to_mesh_file above
     205#domain.set_quantity('elevation',
     206#                     #filename = project.combineddemname + '.pts',
     207#                     filename = project.finedemname + '.pts',
     208                     #filename = project.datadir + 'bathy_land25_orig.pts',
     209#                     use_cache = True,
     210#                     verbose = True)
     211
    173212
    174213# for switching between overall clipping regions - NOT USED
Note: See TracChangeset for help on using the changeset viewer.