source: branches/anuga_1_1/anuga_work/production/west_tas_2008/run_west_tas_smf.py @ 7799

"""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_smf.py
The output sww file is stored in project_smf.outputtimedir 

The scenario is defined by a triangular mesh created from project_smf.polygon,
the elevation data and a tsunami wave generated by s submarine mass failure.

Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
"""

#-------------------------------------------------------------------------------
# Import necessary modules
#-------------------------------------------------------------------------------

# Standard modules
import os
import time
from shutil import copy
from os.path import dirname, basename
from os import mkdir, access, F_OK, sep
import sys

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
from anuga.geospatial_data.geospatial_data import *
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier

# Application specific imports
import project_smf              # Definition of file names and polygons

#-------------------------------------------------------------------------------
# Copy scripts to time stamped output directory and capture screen
# output to file
#-------------------------------------------------------------------------------

# creates copy of code in output dir 
#copy_code_files(project_smf.outputtimedir,__file__,dirname(project_smf.__file__)+sep+ project_smf.__name__+'.py' )
myid = 0
numprocs = 1
#start_screen_catcher(project_smf.outputtimedir, myid, numprocs)
#barrier()

print 'USER:    ', project_smf.user

#-------------------------------------------------------------------------------
# Preparation of topographic data
#
# Convert ASC 2 DEM 2 PTS using source data and store result in source data
#-------------------------------------------------------------------------------

# filenames
onshore_250_dem_name = project_smf.onshore_250_dem_name
meshname = project_smf.meshname+'.msh'

# creates DEM from asc data
convert_dem_from_ascii2netcdf(onshore_250_dem_name, use_cache=True, verbose=True)

#creates pts file for onshore DEM
dem2pts(onshore_250_dem_name, use_cache=True, verbose=True)

print 'create offshore'
G1 = Geospatial_data(file_name = project_smf.offshore_dem_name1 + '.txt')+\
     Geospatial_data(file_name = project_smf.offshore_dem_name2 + '.txt')+\
     Geospatial_data(file_name = project_smf.offshore_dem_name3 + '.txt')

print 'create onshore'
G2 = Geospatial_data(file_name = project_smf.onshore_250_dem_name + '.pts')

print 'create coastline'
G3 = Geospatial_data(file_name = project_smf.coast_line + '.txt')

print 'add'
G = G1.clip(Geospatial_data(project_smf.polyAll)) + G2 + G3

print 'export points'
G.export_points_file(project_smf.combined_dem_name + '.pts')


#----------------------------------------------------------------------------
# Create the triangular mesh based on overall clipping polygon with a tagged
# boundary and interior regions defined in project_smf.py along with
# resolutions (maximal area of per triangle) for each polygon
#-------------------------------------------------------------------------------

from anuga.pmesh.mesh_interface import create_mesh_from_regions
remainder_res = 250000.
region_res = 50000.
local_res = 500.
interior_regions = [[project_smf.poly_region, region_res],
                    [project_smf.poly_local, local_res]]

from caching import cache
_ = cache(create_mesh_from_regions,
          project_smf.polyAll,
           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
                              'e3': [3], 'e4':[4]},
           'maximum_triangle_area': remainder_res,
           'filename': meshname,
           'interior_regions': interior_regions},
          verbose = True, evaluate=False)
print 'created mesh'

#-------------------------------------------------------------------------------                                 
# Setup computational domain
#-------------------------------------------------------------------------------                                 
domain = Domain(meshname, use_cache = True, verbose = True)

print 'Number of triangles = ', len(domain)
print 'The extent is ', domain.get_extent()
print domain.statistics()

domain.set_name(project_smf.basename)
domain.set_datadir(project_smf.outputtimedir)
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)
domain.set_store_vertices_uniquely(False)


#-------------------------------------------------------------------------------                                 
# Setup initial conditions
#-------------------------------------------------------------------------------

tide = 0.0
domain.set_quantity('stage', tide)
domain.set_quantity('friction', 0.0) 
domain.set_quantity('elevation', 
                    filename = project_smf.combined_dem_name + '.pts',
                    use_cache = True,
                    verbose = True,
                    alpha = 0.1
                    )

#-------------------------------------------------------------------------------
# Set up scenario (tsunami_source is a callable object used with set_quantity)
#-------------------------------------------------------------------------------
from smf import slide_tsunami

tsunami_source = slide_tsunami(length=project_smf.length,
                               width=project_smf.width,
                               depth=project_smf.depth,
                               slope=project_smf.slope,
                               thickness=project_smf.thickness, 
                               x0=project_smf.smf_origin[0], 
                               y0=project_smf.smf_origin[1], 
                               alpha=project_smf.alpha, 
                               domain=domain)

#-------------------------------------------------------------------------------                                 
# Setup boundary conditions
#-------------------------------------------------------------------------------
print 'Available boundary tags', domain.get_boundary_tags()

Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])

domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )


#-------------------------------------------------------------------------------                                 
# Evolve system through time
#-------------------------------------------------------------------------------
import time
from Numeric import allclose
from anuga.abstract_2d_finite_volumes.quantity import Quantity

t0 = time.time()

for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'e2')
    stagestep = domain.get_quantity('stage') 

    if allclose(t, 30):
        smf = Quantity(domain)
        smf.set_values(tsunami_source)
        domain.set_quantity('stage', smf + stagestep)
    
print 'That took %.2f seconds' %(time.time()-t0)

print 'finished'