source: production/sydney_2006/run_sydney.py @ 3330

"""Script for running a tsunami inundation scenario for Sydney, NSW, Australia.

Source data such as elevation and boundary data is assumed to be available in
directories specified by project.py
The output sww file is stored in project.outputdir 

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and boundary data obtained from a tsunami simulation done with MOST.

Ole Nielsen, GA - 2005 and Jane Sexton, GA - 2006
"""

tide = 0       #Australian Height Datum (mean sea level)

import os
import time


from pyvolution.shallow_water import Domain, Reflective_boundary, File_boundary,\
     Dirichlet_boundary, Time_boundary, Transmissive_boundary
from pyvolution.data_manager import convert_dem_from_ascii2netcdf,\
     dem2pts
from pyvolution.pmesh2domain import pmesh_to_domain_instance 
from caching import cache
import project
from math import pi, sin

#Data preparation
#Convert ASC 2 DEM 2 PTS using source data and store result in source data
demname = project.demname
meshname = project.meshname+'.msh'

cache(convert_dem_from_ascii2netcdf, demname, {'verbose': True},
      dependencies = [demname + '.asc'],
      verbose = True)
      #evaluate = True)

cache(dem2pts, demname, {'verbose': True},
      dependencies = [demname + '.dem'],      
      verbose = True)

#this allows the user to switch between different clipping polygons
print 'Which total zone are you interested in?'
mytest = int(raw_input('0 = all, 1 = harbour and 2 = botany bay     '))

#Create Triangular Mesh
from pmesh.create_mesh import create_mesh_from_regions

if mytest == 0:
    # for whole region
    south = project.south 
    north = project.north
    west = project.west
    east = project.east

    interior_regions = [[project.harbour_polygon, 25000],
                        [project.botanybay_polygon, 25000]] # maximal area of per triangle
                       
    m = cache(create_mesh_from_regions,
              project.polygonall,
              {'boundary_tags': {'bottom': [0], 'top': [2],
                                 'right': [1], 'left': [3]},
               'resolution': 100000,
               'filename': meshname,           
               'interior_regions': interior_regions},
              evaluate=True,    
              verbose = True)
    #import sys; sys.exit()

if mytest == 1:
    # for harbour region
    south = project.hsouth 
    north = project.hnorth
    west = project.hwest
    east = project.heast

    #interior_regions = [[project.harbour_polygon, 25000]] # maximal area of per triangle
                       
    m = cache(create_mesh_from_regions,
              project.polygon_h,
              {'boundary_tags': {'bottom': [0], 'top': [2],
                                 'right': [1], 'left': [3]},
               'resolution': 50000,
               'filename': meshname},           
     #          'interior_regions': interior_regions},
              evaluate=True, 
              verbose = True)

if mytest == 2:
    # for botany bay region
    south = project.bsouth 
    north = project.bnorth
    west = project.bwest
    east = project.beast

    #interior_regions = [[project.botanybay_polygon, 25000]] # maximal area of per triangle
                       
    m = cache(create_mesh_from_regions,
              project.polygon_bb,
              {'boundary_tags': {'bottom': [0], 'top': [2],
                                 'right': [1], 'left': [3]},
               'resolution': 50000,
               'filename': meshname},
     #          'interior_regions': interior_regions},
              evaluate=True, 
              verbose = True)
    
#Setup domain

domain = cache(pmesh_to_domain_instance, (meshname, Domain),
               dependencies = [meshname],                     
               verbose = True)               


domain.set_name(project.basename)
domain.set_datadir(project.outputdir)
domain.store = True

#domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
domain.quantities_to_be_stored = ['stage']
        
print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()



#Setup Initial Conditions
domain.set_quantity('friction', 0)
domain.set_quantity('stage', tide)
domain.set_quantity('elevation',
                    filename = demname + '.pts',
                    use_cache = True,
                    verbose = True)



#Setup Boundary Conditions
print domain.get_boundary_tags()

Br = Reflective_boundary(domain)
Bt = Transmissive_boundary(domain)
Bd = Dirichlet_boundary([0,0,0])
# 10 min square wave starting at 1 min, 6m high
Bw = Time_boundary(domain=domain,
                   f=lambda t: [(60<t<660)*6, 0, 0])


domain.set_boundary( {'top': Br, 'bottom': Br, 'right': Bw, 'left': Br} )


#Evolve
import time
t0 = time.time()

for t in domain.evolve(yieldstep = 1, finaltime = 3600):
    domain.write_time()
    domain.write_boundary_statistics(tags = 'right') #quantities = 'stage')       

print 'That took %.2f seconds' %(time.time()-t0)