source: anuga_work/production/wollongong_2006/run_kembla.py @ 7855

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

Source data such as elevation and boundary data is assumed to be available in
directories specified by project_kembla.py

The scenario is defined by a triangular mesh created from project_kembla.polygon,
the elevation data and an input wave.

Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis and Jane Sexton, 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, Time_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

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

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

# filenames
dem_name = project_kembla.dem_name
meshname = project_kembla.meshname+'.msh'

### creates DEM from asc data
##convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)
##
###creates pts file for onshore DEM
##dem2pts(dem_name, use_cache=True, verbose=True)
##
##print 'create onshore'
##G = Geospatial_data(file_name = project_kembla.dem_name + '.pts')
##
##print 'export points'
##G.export_points_file(project_kembla.combined_dem_name + '.pts')
#G.export_points_file(project_kembla.combined_dem_name + '.xya')

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

from anuga.pmesh.mesh_interface import create_mesh_from_regions
remainder_res = 100000
bay_res = 5000
interior_regions = [[project_kembla.poly_bay, bay_res]]

from caching import cache
_ = cache(create_mesh_from_regions,
          project_kembla.polyAll,
           {'boundary_tags': {'top': [0], 'right': [1], 'bottom': [2],
                              'left': [3]},
           '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_kembla.basename)
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)

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

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

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

Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,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': Bd,'right': Bw, 'bottom': Bd, 'left': Bd} )


#-------------------------------------------------------------------------------                                 
# Evolve system through time
#-------------------------------------------------------------------------------
import time
t0 = time.time()
from Numeric import allclose
from anuga.abstract_2d_finite_volumes.quantity import Quantity

for t in domain.evolve(yieldstep = 10, finaltime = 1000): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'right')
    stagestep = domain.get_quantity('stage') 
    
print 'That took %.2f seconds' %(time.time()-t0)

print 'finished'