Changeset 2640


Ignore:
Timestamp:
Mar 31, 2006, 1:29:48 PM (18 years ago)
Author:
sexton
Message:

Update Sydney example post Benfield visit (30 March 2006) to mimic changes from island.py

File:
1 edited

Legend:

Unmodified
Added
Removed
  • production/sydney_2006/run_sydney_smf.py

    r2537 r2640  
    2020
    2121# Related major packages
    22 from pyvolution.shallow_water import Domain, Reflective_boundary
     22from pyvolution.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary
    2323from pyvolution.data_manager import convert_dem_from_ascii2netcdf, dem2pts
     24#from pyvolution.data_manager_old import convert_dem_from_ascii2netcdf, dem2pts
    2425from pyvolution.combine_pts import combine_rectangular_points_files
    2526from pyvolution.pmesh2domain import pmesh_to_domain_instance
     27from pyvolution.quantity import Quantity
     28from Numeric import allclose
    2629
    2730# Application specific imports
     
    5760        verbose=True)
    5861
    59 
    6062# combining the coarse and fine data
    6163combine_rectangular_points_files(project.finedemname + '.pts',
     
    6466
    6567#from pmesh.create_mesh import create_mesh_from_regions
    66 # new interface
     68#new interface
    6769from pmesh.mesh_interface import create_mesh_from_regions
    6870
    69 # original
    70 interior_res = 5000
    71 interior_regions = [[project.harbour_polygon_2, interior_res],
    72                     [project.botanybay_polygon_2, interior_res]]
     71# original issue to Benfield
     72#interior_res = 5000
     73#interior_regions = [[project.harbour_polygon_2, interior_res],
     74#                    [project.botanybay_polygon_2, interior_res]]
     75
     76# used for finer mesh
     77interior_res1 = 5000
     78interior_res2 = 315
     79interior_regions = [[project.newpoly1, interior_res1],
     80                    [project.south1, interior_res1],
     81                    [project.finepolymanly, interior_res2],
     82                    [project.finepolyquay, interior_res2]]
     83
     84print 'number of interior regions', len(interior_regions)
    7385
    7486#FIXME: Fix caching of this one once the mesh_interface is ready
    7587from caching import cache
    7688
     89# original + refined region
    7790_ = cache(create_mesh_from_regions,
    7891          project.diffpolygonall,
     
    8699          verbose = True)
    87100
    88 
    89 #-----------------------------------------------------------------------------
     101#-------------------------------------------------------------------------------                                 
    90102# Setup computational domain
    91 #-----------------------------------------------------------------------------
     103#-------------------------------------------------------------------------------                               
    92104
    93105domain = pmesh_to_domain_instance(meshname, Domain,
     
    97109print 'Number of triangles = ', len(domain)
    98110print 'The extent is ', domain.get_extent()
     111print domain.statistics()
    99112
    100113domain.set_name(project.basename)
     
    119132
    120133
    121 #------------------------------------------------------------------------------                                 
     134#-------------------------------------------------------------------------------                                 
    122135# Setup initial conditions
    123136#-------------------------------------------------------------------------------
    124137
    125 domain.set_quantity('stage', tsunami_source)
    126 domain.set_quantity('friction', 0.0)
     138# apply slump after 30 mins, initialise to water level of tide = 0
     139domain.set_quantity('stage', 0.0)
     140domain.set_quantity('friction', 0.03)
    127141domain.set_quantity('elevation',
    128142                    filename = project.combineddemname + '.pts',
     
    138152
    139153Br = Reflective_boundary(domain)
    140 domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
    141                       'right2': Br, 'top': Br, 'left1': Br,
    142                       'left2': Br, 'left3': Br} )
     154Bd = Dirichlet_boundary([0, 0, 0])
     155
     156# original + refined regions
     157#domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
     158#                      'right2': Br, 'top': Br, 'left1': Br,
     159#                      'left2': Br, 'left3': Br} )
     160
     161# enforce Dirichlet BC - from 30/03/06 Benfield visit
     162domain.set_boundary( {'bottom': Bd, 'right1': Bd, 'right0': Bd,
     163                      'right2': Bd, 'top': Bd, 'left1': Bd,
     164                      'left2': Bd, 'left3': Bd} )
     165
     166# increasingly finer interior regions
     167#domain.set_boundary( {'bottom': Br, 'right': Br, 'left': Br, 'top': Br} )
    143168
    144169
     
    149174import time
    150175t0 = time.time()
     176thisfile = project.integraltimeseries+'.csv'
     177fid = open(thisfile, 'w')
    151178
    152 for t in domain.evolve(yieldstep = 120, finaltime = 18000): 
     179for t in domain.evolve(yieldstep = 120, finaltime = 18000):
    153180    domain.write_time()
    154     domain.write_boundary_statistics(tags = 'bottom')     
     181    domain.write_boundary_statistics(tags = 'bottom')
     182
     183    # calculate integral
     184    thisstagestep = domain.get_quantity('stage')
     185    s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral())
     186    fid.write(s)
     187
     188    # add slump after 30 mins
     189    if allclose(t, 30*60):
     190        slump = Quantity(domain)
     191        slump.set_values(tsunami_source)
     192        domain.set_quantity('stage', slump + thisstagestep)
     193   
    155194   
    156195print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset for help on using the changeset viewer.