source: production/sydney_2006/run_sydney_smf.py @ 2292

"""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 Adrian Hitchman 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 pyvolution.combine_pts import combine_rectangular_points_files 
from caching import cache
import project
from math import pi, sin
from smf import slump_tsunami, slide_tsunami, Double_gaussian
from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 

# Data preparation
# Convert ASC 2 DEM 2 PTS using source data and store result in source data
# Do for coarse and fine data
# Fine pts file to be clipped to area of interest
coarsedemname = project.coarsedemname
finedemname = project.finedemname
meshname = project.meshname+'.msh'

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

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

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

# clipping the fine data
cache(dem2pts, finedemname, {'verbose': True,
      'easting_min': project.eastingmin,
      'easting_max': project.eastingmax,
      'northing_min': project.northingmin,
      'northing_max': project.northingmax},
      dependencies = [finedemname + '.dem'],
      #evaluate = True,
      verbose = True)

# combining the coarse and fine data
combine_rectangular_points_files(project.finedemname + '.pts',
                                 project.coarsedemname + '.pts',
                                 project.combineddemname + '.pts')
                                 
# 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     '))

mytest = 0

# Create Triangular Mesh
# Overall clipping polygon and interior regions defined in project.py
from pmesh.create_mesh import create_mesh_from_regions

if mytest == 0:
    # for whole region
    interior_res = 2000 # run at 2000 for final one
    interior_regions = [[project.harbour_polygon_2, interior_res],
                        [project.botanybay_polygon_2, interior_res]] # maximal area of per triangle

    m = cache(create_mesh_from_regions,
              project.diffpolygonall,
              {'boundary_tags': {'bottom': [0], 
                                 'right1': [1], 'right0': [2],
                                 'right2': [3], 'top': [4], 'left1': [5],
                                 'left2': [6], 'left3': [7]},
               'resolution': 100000,
               'filename': meshname,           
               'interior_regions': interior_regions},
              #evaluate=True,    
              verbose = True)
    #import sys; sys.exit()
    
    #Add elevation data to the mesh - this is in place of set_quantity
    mesh_elevname = 'elevation.msh'   
    cache(fit_to_mesh_file,(meshname,
                     #project.finedemname + '.pts',
                     #project.coarsedemname + '.pts',
                     project.combineddemname + '.pts',
                     mesh_elevname,
                     DEFAULT_ALPHA),
            {'verbose': True, 'expand_search': True, 'precrop': True},
          dependencies = [meshname],
          #evaluate = True,   
          verbose = False)
    meshname = mesh_elevname
    
# Setup domain
domain = cache(pmesh_to_domain_instance, (meshname, Domain),
               dependencies = [meshname],                     
               verbose = True)               
        
print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()

domain.set_name(project.basename)
domain.set_datadir(project.outputdir)
domain.store = True
#domain.quantities_to_be_stored = ['stage']
domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']

#print 'Do you want to run the slump scenario?'
#q2 = int(raw_input('0 for yes, 1 for no     '))

q2 = 0

if q2 == 0:
    # slump parameters
    # as advised by Adrian, this can be as much as 15000-30000 metres
    # len = 4500.0
    len = 15000.0
    thk = 670 #thk = 760.0
    wid = 15000.0 # for slump scenario, wid=len
    #dep = 400 #1000.0 #d = 400 doesn't make sense with thk = 670!
    dep = 1000.0
    rad = 3330
    dp = 0.23
    th = 6
    alpha = 0.0
    x0 = project.x0
    y0 = project.y0
    # slide parameters
    # len = 10000.0
    # dep = 1200.0
    # thk = 200.0
    # th = 12
    # Slump or Slide Initial conditions
    t = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
                      radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
    #t = slide_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
    #                  x0=x0, y0=y0, alpha=alpha)
    domain.set_quantity('stage', t, verbose = True)

if q2 == 1:
    domain.set_quantity('stage', tide)


# These quantities are the same regardless of the initial condition
domain.set_quantity('friction', 0)

    
# 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: [(6<t<606)*6, 0, 0])

# for slump scenario
if q2 == 0:
    domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
                                 'right2': Br, 'top': Br, 'left1': Br,
                                 'left2': Br, 'left3': Br} )

# for square wave scenario
if q2 == 1:
    domain.set_boundary( {'top': Br, 'bottom': Br, 'right': Bw, 'left': Br} )

# Evolve
import time
t0 = time.time()

for t in domain.evolve(yieldstep = 120, finaltime = 10800): #120 10800 = 3hrs
    domain.write_time()
    domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage')     
    
print 'That took %.2f seconds' %(time.time()-t0)

#############################################################
# testing system on a smaller mesh
    #m = create_mesh_from_regions(project.harbour_polygon,
    #                             boundary_tags={'bottom': [0], 'top': [2],
    #                                            'right': [1], 'left': [3]},
    #                             resolution=100000, filename=meshname)
    
# This section replaced with fit_to_mesh_file above
#domain.set_quantity('elevation',
#                     #filename = project.combineddemname + '.pts',
#                     filename = project.finedemname + '.pts',
                     #filename = project.datadir + 'bathy_land25_orig.pts',
#                     use_cache = True,
#                     verbose = True)


# for switching between overall clipping regions - NOT USED

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
    # meshname = project.meshname + 'harbour.msh'                   
    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
    # meshname = project.meshname + 'bay.msh'                   
    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)