source: anuga_work/production/australia_ph2/sydney/comparisons/run_sydney_250m.py @ 6295

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

The scenario is defined by a triangular mesh created from project_250m.polygon,
the elevation data is compiled into a pts file through build_perth.py
and a simulated tsunami is generated through an sts file from build_boundary.py.

Input: sts file (build_boundary.py for respective event)
       pts file (build_perth.py)
       information from project file
Outputs: sww file stored in project_250m.output_run_time_dir 
The export_results_all.py and get_timeseries.py is reliant
on the outputs of this script

Ole Nielsen and Duncan Gray, GA - 2005, Jane Sexton, Nick Bartzis, GA - 2006
Ole Nielsen, Jane Sexton and Kristy Van Putten - 2008
"""

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

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

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water.shallow_water_domain import Transmissive_stage_zero_momentum_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import File_boundary
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Field_boundary
from anuga.shallow_water import Time_boundary
from Numeric import allclose
from anuga.shallow_water.data_manager import export_grid, create_sts_boundary
from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
from anuga_parallel.parallel_abstraction import get_processor_name
from anuga.caching import myhash
from anuga.damage_modelling.inundation_damage import add_depth_and_momentum2csv, inundation_damage
from anuga.fit_interpolate.benchmark_least_squares import mem_usage 
from anuga.utilities.polygon import read_polygon, plot_polygons, polygon_area, is_inside_polygon
from anuga.geospatial_data.geospatial_data import find_optimal_smoothing_parameter
from polygon import Polygon_function
    
# Application specific imports
import project_250m  # Definition of file names and polygons
numprocs = 1
myid = 0

def run_model(**kwargs):
    
    #------------------------------------------------------------------------------
    # Copy scripts to time stamped output directory and capture screen
    # output to file
    #------------------------------------------------------------------------------
    print "Processor Name:",get_processor_name()

    #copy script must be before screen_catcher

    print 'output_dir',kwargs['output_dir']
    
    copy_code_files(kwargs['output_dir'],__file__, 
             dirname(project_250m.__file__)+sep+ project_250m.__name__+'.py' )

    store_parameters(**kwargs)

    start_screen_catcher(kwargs['output_dir'], myid, numprocs)

    print "Processor Name:",get_processor_name()
    
    #-----------------------------------------------------------------------
    # Domain definitions
    #-----------------------------------------------------------------------
##
##    # Read in boundary from ordered sts file
##    urs_bounding_polygon=create_sts_boundary(os.path.join(project_250m.boundaries_dir_event,project_250m.scenario_name))
##
##    # Reading the landward defined points, this incorporates the original clipping
##    # polygon minus the 100m contour
##    landward_bounding_polygon = read_polygon(project_250m.landward_dir)
##
##    # Combine sts polyline with landward points
##    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
##    
##    # counting segments
##    N = len(urs_bounding_polygon)-1
##    print 'N'+str(N)
##    # boundary tags refer to project_250m.landward 4 points equals 5 segments start at N
##    boundary_tags={'back': [N+1,N+2,N+3,N+4,N+5,N+6,N+7,N+8,N+9,N+10,N+11,N+12,N+13,N+14,N+15,N+16,N+17,N+18,N+19], 'side': [N,N+20], 'ocean': range(N)}
    if project_250m.area =='large':
        land = [2,3,4,]
        sea = [0]
        side1 = [1]
        side2 = [5]
        bounding_polygon = project_250m.poly_all
        boundary_tags={'back': land, 'ocean': sea, 'side1':side1,'side2':side2}
    elif project_250m.area =='small':
        land = [2,3,4,]
        sea = [0]
        side1 = [1]
        side2 = [5]
        bounding_polygon = project_250m.poly_all
        boundary_tags={'back': land, 'ocean': sea, 'side1':side1,'side2':side2}
    elif project_250m.area =='medium':
        land = [2,3,4]
        sea = [0]
        side1 = [1]
        side2 = [5]
        bounding_polygon = project_250m.poly_all
        boundary_tags={'back': land, 'ocean': sea, 'side1':side1,'side2':side2}
    elif project_250m.area =='vlarge':
        land = [9,10,11,12,13,14]
        sea = [0,1,2,3,4,5,6,7]
        side1 = [8]
        side2 = [15]
        bounding_polygon = project_250m.poly_all
        boundary_tags={'back': land, 'ocean': sea, 'side1':side1,'side2':side2} 
    else:
        print 'area not defined therefore bondary_tags not defined'
    #--------------------------------------------------------------------------
    # Create the triangular mesh based on overall clipping polygon with a tagged
    # boundary and interior regions defined in project_250m.py along with
    # resolutions (maximal area of per triangle) for each polygon
    #--------------------------------------------------------------------------

    # IMPORTANT don't cache create_mesh_from_region and Domain(mesh....) as it
    # causes problems with the ability to cache set quantity which takes alot of times
        
    print 'start create mesh from regions'

    create_mesh_from_regions(bounding_polygon,
                         boundary_tags=boundary_tags,
                         maximum_triangle_area=project_250m.res_poly_all,
                         interior_regions=project_250m.interior_regions,
                         filename=project_250m.meshes_dir_name,
                         use_cache=False,
                         verbose=True)
   
    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
    print 'Setup computational domain'

    domain = Domain(project_250m.meshes_dir_name, use_cache=False, verbose=True)
    print 'memory usage before del domain',mem_usage()
       
    print domain.statistics()
    print 'triangles',len(domain)
    
    kwargs['act_num_trigs']=len(domain)


    #-------------------------------------------------------------------------
    # Setup initial conditions
    #-------------------------------------------------------------------------
    print 'Setup initial conditions'

    # sets the initial stage in the offcoast region only
    IC = Polygon_function( [(project_250m.poly_mainland, 0),(project_250m.island1,0),(project_250m.island2,0)], default = kwargs['tide'],
                             geo_reference = domain.geo_reference)
    domain.set_quantity('stage', IC)
    #domain.set_quantity('stage',kwargs['tide'] )
    domain.set_quantity('friction', kwargs['friction']) 
    
    print 'Start Set quantity',kwargs['elevation_file']

    domain.set_quantity('elevation', 
                        filename = kwargs['elevation_file'],
                        use_cache = False,
                        verbose = True,
                        alpha = kwargs['alpha'])
    print 'Finished Set quantity'

##   #------------------------------------------------------
##    # Distribute domain to implement parallelism !!!
##    #------------------------------------------------------
##
##    if numprocs > 1:
##        domain=distribute(domain)

    #------------------------------------------------------
    # Set domain parameters
    #------------------------------------------------------
    print 'domain id', id(domain)
    domain.set_name(kwargs['scenario_name'])
    domain.set_datadir(kwargs['output_dir'])
    domain.set_default_order(2)                 # Apply second order scheme
    domain.set_minimum_storable_height(0.01)    # Don't store anything less than 1cm
    domain.set_store_vertices_uniquely(False)
    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    domain.tight_slope_limiters = 1
    print 'domain id', id(domain)

    #-------------------------------------------------------------------------
    # Setup boundary conditions 
    #-------------------------------------------------------------------------
    print 'Available boundary tags', domain.get_boundary_tags()
    print 'domain id', id(domain)
   
    boundary_urs_out=project_250m.boundaries_dir_event + sep + project_250m.scenario_name

    Br = Reflective_boundary(domain)
    Bd = Dirichlet_boundary([kwargs['tide'],0,0])
    Bt = Transmissive_stage_zero_momentum_boundary(domain)
    Rb = Time_boundary (domain, f=lambda t: [(sin(t*pi/600) * exp(-(t/3600)**2)), 0.0, 0.0]) 
##    print 'Available boundary tags', domain.get_boundary_tags()
##    Bf = Field_boundary(boundary_urs_out+'.sts',  # Change from file_boundary
##                   domain, mean_stage= project_250m.tide,
##                   time_thinning=1,
##                   default_boundary=Bd,
##                   use_cache=False,
##                   verbose = True,
##                   boundary_polygon=bounding_polygon)

    domain.set_boundary({'back': Br,
                         'side1': Bt,
                         'side2':Bt,
                         'ocean':Rb}) 

    kwargs['input_start_time']=domain.starttime

    print'finish set boundary'

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

    for t in domain.evolve(yieldstep = project_250m.yieldstep, finaltime = kwargs['finaltime']
                       ,skip_initial_step = False): 
        domain.write_time()
        domain.write_boundary_statistics(tags = 'ocean')

    # these outputs should be checked with the resultant inundation map
    x, y = domain.get_maximum_inundation_location()
    q = domain.get_maximum_inundation_elevation()
    print 'Maximum runup observed at (%.2f, %.2f) with elevation %.2f' %(x,y,q)

    print 'Simulation took %.2f seconds' %(time.time()-t0)

    #kwargs 'completed' must be added to write the final parameters to file
    kwargs['completed']=str(time.time()-t0)
     
    store_parameters(**kwargs)
     
    print 'memory usage before del domain1',mem_usage()
    
    
#-------------------------------------------------------------
if __name__ == "__main__":
    
    kwargs={}
    kwargs['finaltime']=project_250m.finaltime
    kwargs['output_dir']=project_250m.output_run_time_dir
    kwargs['elevation_file']=project_250m.combined_dir_name+'.pts'
    kwargs['scenario_name']=project_250m.scenario_name
    kwargs['tide']=project_250m.tide
    kwargs['alpha'] = project_250m.alpha
    kwargs['friction']=project_250m.friction
     
    run_model(**kwargs)