source: anuga_work/production/MOST_example/run_pt_hedland.py @ 4631

"""Script for running a tsunami inundation scenario for Onslow, WA, Australia.

 Boundary data is assumed to be in this directory.
Move the SU-AU_clip_blahblah.nc files into this directory, from teh pt_hedland
data directory.

The file cripple_pts_file creates the elevation file
pt_hedland_combined_elevation_31204.pts, from
pt_hedland_combined_elevation.pts - which is created by the other
run_pt_hedland script.
 
directories specified by project.py
The output sww file is stored in project.outputtimedir 

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and a simulated submarine landslide.

Use shallow_water.most2nc to create the SU-AU_clip_e.nc file.

Ole Nielsen, Duncan Gray and Nick Bartzis, GA - 2005,2006
"""
#-------------------------------------------------------------------------------# Import necessary modules
#-------------------------------------------------------------------------------

# Standard modules
from os import sep
from os.path import dirname, basename
import time

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, \
                            Dirichlet_boundary, Time_boundary, File_boundary
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, \
     dem2pts, ferret2sww
from shutil import copy
from os import mkdir, access, F_OK
from anuga.geospatial_data.geospatial_data import *
import sys
from anuga.abstract_2d_finite_volumes.util import Screen_Catcher

from anuga.fit_interpolate.fit import fit_to_mesh_file

# Application specific imports
import project                 # 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 if dir doesn't exist
if access(project.outputtimedir,F_OK) == 0 :
    mkdir (project.outputtimedir)
copy (dirname(project.__file__) +sep+ project.__name__+'.py', project.outputtimedir + project.__name__+'.py')
copy (__file__, project.outputtimedir + basename(__file__))
print 'project.outputtimedir',project.outputtimedir

# normal screen output is stored in 
#screen_output_name = project.outputtimedir + "screen_output.txt"
#screen_error_name = project.outputtimedir + "screen_error.txt"

# used to catch screen output to file
#sys.stdout = Screen_Catcher(screen_output_name)
#sys.stderr = Screen_Catcher(screen_error_name)
print 'USER:    ', project.user

#-------------------------------------------------------------------------------
# Preparation of topographic data
#
# 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
#-------------------------------------------------------------------------------

# filenames
meshname = project.meshname + '.tsh'
mesh_elevname = project.mesh_elevname  + '.tsh'
source_dir = project.boundarydir


#-------------------------------------------------------------------------------                                 
# 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

region_res = 5000000
coast_res = 500000
pt_hedland_res = 500000
interior_regions = [[project.poly_pt_hedland, pt_hedland_res],
                    [project.poly_region, region_res]]

print 'number of interior regions', len(interior_regions)

print 'start create mesh from regions'
from caching import cache
_ = cache(create_mesh_from_regions,
          project.polyAll,
          {'boundary_tags': {'topright': [0], 'topleft': [1],
                             'left': [2], 'bottom0': [3],
                             'bottom1': [4], 'bottom2': [5],
                             'bottom3': [6], 'right': [7]},
           'maximum_triangle_area': 5000000,
           'filename': meshname,           
           'interior_regions': interior_regions},
          verbose = True
          #, evaluate=True
          )

cache(fit_to_mesh_file,(meshname,
                 'pt_hedland_combined_elevation_31204' + '.pts',
                 mesh_elevname),
      {'verbose': True}
      #,evaluate = True      
      ,verbose = False
      )
#-------------------------------------------------------------------------------                                 
# Setup computational domain
#-------------------------------------------------------------------------------                                 
domain = Domain(mesh_elevname, use_cache = False, verbose = True)

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

domain.set_name(project.basename)
domain.set_datadir(project.outputtimedir)
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])

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

tide = 0.
#high
#tide = 3.6
#low
#tide = -3.9

domain.set_quantity('stage', tide)
domain.set_quantity('friction', 0.0) 
print 'hi and file',project.combined_dem_name + '.pts'

#domain.set_quantity('elevation', 
#                    filename = project.combined_dem_name + '.pts',
#                    use_cache = True,
#                    verbose = True,
#                    alpha = 0.1
 #                   )

#-------------------------------------------------------------------------------                                 
# Setup boundary conditions (all reflective)
#-------------------------------------------------------------------------------
print 'start ferret2sww'
# skipped as results in file SU-AU_clipped is correct for all WA


south = project.south 
north = project.north
west = project.west
east = project.east

#note only need to do when an SWW file for the MOST boundary doesn't exist
cache(ferret2sww,
      (project.boundary_basename,
       project.boundary_basename+'_'+project.basename), 
      {'verbose': True,
       'minlat': south,
       'maxlat': north,
       'minlon': west,
       'maxlon': east,
#       'origin': project.mesh_origin,
       'origin': domain.geo_reference.get_origin(),
       'mean_stage': tide,
       'zscale': 10,                 #Enhance tsunami
       'fail_on_NaN': False,
       'inverted_bathymetry': True},
       #evaluate = True,
       verbose = True,
      dependencies = source_dir + project.boundary_basename + '.sww')

print 'Available boundary tags', domain.get_boundary_tags()

Bf = File_boundary(project.boundary_basename+'_'+project.basename + '.sww', 
                    domain, verbose = True)
Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])
domain.set_boundary( {'topright': Bf,'topleft': Bf, 'left':  Bd, 'bottom0': Bd,
                      'bottom1': Bd, 'bottom2': Bd, 'bottom3': Bd, 
                        'right': Bd})

#-------------------------------------------------------------------------------                                 
# Evolve system through time
#-------------------------------------------------------------------------------
import time
t0 = time.time()

for t in domain.evolve(yieldstep = 240, finaltime = 10800): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'topright')     

for t in domain.evolve(yieldstep = 120, finaltime = 16200
                       ,skip_initial_step = True): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'topright')     

for t in domain.evolve(yieldstep = 60, finaltime = 21600
                       ,skip_initial_step = True): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'topright')     
    
for t in domain.evolve(yieldstep = 120, finaltime = 27000
                       ,skip_initial_step = True): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'topright')     

for t in domain.evolve(yieldstep = 240, finaltime = 36000
                       ,skip_initial_step = True): 
    domain.write_time()
    domain.write_boundary_statistics(tags = 'topright')   
  
print 'That took %.2f seconds' %(time.time()-t0)

print 'finished'