source: anuga_work/production/patong/run_patong.py @ 6133

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

The scenario is defined by a triangular mesh created from project.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)
       pts file (build_perth.py)
       information from project file
Outputs: sww file stored in project.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
import time
import sys

# Related major packages
from anuga.shallow_water import Domain
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 Numeric import allclose
from anuga.shallow_water.data_manager import export_grid, create_sts_boundary, csv2building_polygons
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.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  # 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
    #------------------------------------------------------------------------------

    #copy script must be before screen_catcher

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

    store_parameters(**kwargs)

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

    
    #-----------------------------------------------------------------------
    # Domain definitions
    #-----------------------------------------------------------------------

    # Read in boundary from ordered sts file
    urs_bounding_polygon=create_sts_boundary(os.path.join(project.boundaries_dir,project.scenario_name))

    # Reading the landward defined points, this incorporates the original clipping
    # polygon minus the 100m contour
    landward_bounding_polygon = read_polygon(project.landward_dir)

    # Combine sts polyline with landward points
    bounding_polygon = urs_bounding_polygon + landward_bounding_polygon
    
    # counting segments
    N = len(urs_bounding_polygon)-1

    # boundary tags refer to project.landward 4 points equals 5 segments start at N
    boundary_tags={'back': [N+1,N+2,N+3], 'side': [N,N+4], 'ocean': range(N)}

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

    # 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.res_poly_all,
                             interior_regions=project.interior_regions,
                             filename=project.meshes_dir_name,
                             use_cache=True,
                             verbose=True)
   
    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
    print 'Setup computational domain'

    domain = Domain(project.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.poly_mainland, 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'])

    # Add buildings from file
    building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)

    L = []
    for key in building_polygons:
        poly = building_polygons[key]
        elev = building_heights[key]
        L.append((poly, elev))
        
    #Q1 = domain.get_quantity('elevation')
    #Q2 = Quantity(domain) # Temporary quantity for buildings
    #Q2.set_values(Polygon_function(L, default=0.0))
    domain.add_quantity('elevation', Polygon_function(L, default=0.0))
    
    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.boundaries_dir + sep + project.scenario_name

    Br = Reflective_boundary(domain)
    Bd = Dirichlet_boundary([kwargs['tide'],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.tide,
                   time_thinning=1,
                   default_boundary=Bd,
                   use_cache=True,
                   verbose = True,
                   boundary_polygon=bounding_polygon)

    domain.set_boundary({'back': Br,
                         'side': Bd,
                         'ocean': Bf}) 

    kwargs['input_start_time']=domain.starttime

    print'finish set boundary'

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

    for t in domain.evolve(yieldstep = project.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)
     
    
    
#-------------------------------------------------------------
if __name__ == "__main__":
    
    kwargs={}
    kwargs['file_name']=project.dir_comment
    kwargs['finaltime']=project.finaltime
    kwargs['output_dir']=project.output_run_time_dir
    kwargs['elevation_file']=project.combined_dir_name+'.pts'
    kwargs['scenario_name']=project.scenario_name
    kwargs['tide']=project.tide
    kwargs['alpha'] = project.alpha
    kwargs['friction']=project.friction
    #kwargs['num_cpu']=numprocs
    #kwargs['host']=project.host
    #kwargs['starttime']=project.starttime
    #kwargs['yieldstep']=project.yieldstep
    #kwargs['user']=project.user
    #kwargs['time_thinning'] = project.time_thinning
      
    run_model(**kwargs)