source: production/sydney_2006/run_sydney_smf.py @ 2350

"""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')
                                 
# Create Triangular Mesh
# Overall clipping polygon and interior regions defined in project.py
from pmesh.create_mesh import create_mesh_from_regions

# for whole region
interior_res = 5000 # maximal area of per triangle
interior_regions = [[project.harbour_polygon_2, interior_res],
                    [project.botanybay_polygon_2, interior_res]] 

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)

#Add elevation data to the mesh - this is in place of set_quantity
#mesh_elevname = 'elevation.msh'   
#cache(fit_to_mesh_file,(meshname, 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)

# This section replaced with fit_to_mesh_file above
domain.set_quantity('elevation',
                     filename = project.combineddemname + '.pts',
#                     filename = project.finedemname + '.pts',
                     use_cache = True,
                     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', 'xmomentum', 'ymomentum']

# slump parameters, wid=len
len = 30000.0
dep = 400.0
#thk = 120.0 for scenario 1 and 176 for scenario 2 rang 0.2-0.44 of depth
thk = 0.44*dep 
rad = 3330
dp = 0.23
th = 6
alpha = 0.0
x0 = project.x0 # slump origin
y0 = project.y0
# Setup Initial conditions
t = slump_tsunami(length=len, depth=dep, slope=th, thickness=thk, \
                  radius=rad, dphi=dp, x0=x0, y0=y0, alpha=alpha)
domain.set_quantity('stage', t)
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])

domain.set_boundary( {'bottom': Br, 'right1': Br, 'right0': Br,
                      'right2': Br, 'top': Br, 'left1': Br,
                      'left2': Br, 'left3': Br} )

# Evolve
import time
t0 = time.time()

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