source: anuga_work/production/shark_bay_2007/run_shark_bay.py @ 4689

"""Script for running tsunami inundation scenario for Dampier, WA, Australia.

Source data such as elevation and boundary data is assumed to be available in
directories specified by project_urs.py
The output sww file is stored in project_urs.output_run_time_dir 

The scenario is defined by a triangular mesh created from project_urs.polygon,
the elevation data and a simulated tsunami generated with URS code.

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

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

# Standard modules
from os import sep,umask
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 Reflective_boundary
from anuga.shallow_water import Field_boundary
from Numeric import allclose

from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water.data_manager import start_screen_catcher
from anuga.shallow_water.data_manager import copy_code_files
from anuga.shallow_water.data_manager import store_parameters
from anuga_parallel.parallel_api import distribute, numprocs, myid, barrier
from anuga_parallel.parallel_abstraction import get_processor_name
from anuga.caching import myhash

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

def run_model(**kwargs):
    
    tide = kwargs['tide']
    alpha = kwargs['alpha']
    friction = kwargs['friction']
    time_thinning = kwargs['time_thinning']
    scenario_name = kwargs['aa_scenario_name']
    
    #------------------------------------------------------------------------------
    # Copy scripts to time stamped output directory and capture screen
    # output to file
    #------------------------------------------------------------------------------

    #copy script must be before screen_catcher
    print 'tide',tide
    kwargs['est_num_trigs']=project.trigs_min
    kwargs['num_cpu']=numprocs
    kwargs['host']=project.host
    kwargs['res_factor']=project.res_factor
    kwargs['starttime']=project.starttime
    kwargs['yieldstep']=project.yieldstep
    kwargs['finaltime']=project.finaltime
    
    kwargs['output_dir']=project.output_run_time_dir
    kwargs['bathy_file']=project.combined_dir_name + '.pts'
#    kwargs['bathy_file']=project.combined_small_dir_name + '.pts'
    kwargs['boundary_file']=project.boundary_name + '.sww'
#    kwargs['Completed']=''

    print 'output_dir',kwargs['output_dir']
    if myid == 0:
        copy_code_files(kwargs['output_dir'],__file__, 
                        dirname(project.__file__)+sep+project.__name__+'.py' )

        store_parameters(**kwargs)

    barrier()

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

    print 'Processor Name:', get_processor_name()

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

    # FIXME (Ole): Move to build_shark_bay.py
    if myid == 0:
    
        print 'start create mesh from regions'
        
        create_mesh_from_regions(project.bounding_polygon,
                                 boundary_tags=project.boundary_tags,
                                 maximum_triangle_area=project.res_bounding_polygon,
                                 interior_regions=project.interior_regions,
                                 filename=project.mesh_name+'.msh',
                                 use_cache=False,
                                 verbose=True)

    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
    print 'Setup computational domain'

    #domain = cache(Domain, (mesh_name), {'use_cache':True, 'verbose':True}, verbose=True)
    #above don't work
    domain = Domain(project.mesh_name+'.msh',
                    use_cache=False, verbose=True)
     
    print domain.statistics()
    print 'triangles',len(domain)
    
    kwargs['act_num_trigs']=len(domain)

    #-------------------------------------------------------------------------
    # Setup initial conditions
    #-------------------------------------------------------------------------
    if myid == 0:

        IC = project.tide
        #if project.tide == 0.0:
        #    print 'Tide is MSL'
        #    IC = 0.0
        #else:    
        #    print 'Create onshore polygon'
        #    from polygon import Polygon_function
        #    #following sets the stage/water to be offcoast only
        #    IC = Polygon_function([(project.onshore_polygon, -1.0)],
        #                          default = tide,
        #                          geo_reference = domain.geo_reference)
            

        print 'Set initial conditions'
        domain.set_quantity('stage', IC)
        domain.set_quantity('friction', friction)
        print 'Set elevation'
        domain.set_quantity('elevation', 
                            filename = kwargs['bathy_file'],
                            use_cache = True,
                            verbose = True,
                            alpha = alpha)
        

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

    #------------------------------------------------------
    # Set domain parameters
    #------------------------------------------------------
    print 'domain id', id(domain)
    domain.set_name(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.set_maximum_allowed_speed(0.1) # Allow a little runoff (0.1 is OK)
    print 'domain id', id(domain)
    domain.beta_h = 0
    #domain.tight_slope_limiters = 1
    

    #-------------------------------------------------------------------------
    # Setup boundary conditions 
    #-------------------------------------------------------------------------
    print 'Available boundary tags', domain.get_boundary_tags()
    print 'domain id', id(domain)
    #print 'Reading Boundary file',project.boundaries_dir_namea + '.sww'

    Bf = Field_boundary(kwargs['boundary_file'],
                        domain,
                        time_thinning=time_thinning,
                        mean_stage=tide,
                        momentum_scale=project.momentum_scale,                        
                        use_cache=False,
                        verbose=True)
                    
    kwargs['input_start_time']=domain.starttime

    print 'finished reading boundary file'

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

    print'set_boundary'

    #domain.set_boundary({'back': Br,
    #                     'side': Bd,
    #                     'ocean': Bf})
    domain.set_boundary({'tide': Bd,
    #                     'ocean': Dirichlet_boundary([3,0,0])})
                         'ocean': Bf})     
    print'finish set boundary'

    #----------------------------------------------------------------------------
    # Evolve system through time
    #----------------------------------------------------------------------------

    t0 = time.time()

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


    # Reset stage in study area
    print 'Reset initial condition in study area'
    domain.set_quantity('stage', project.tide,
                        polygon=project.poly_tsunami_approach_area)
    

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

    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 'That took %.2f seconds' %(time.time()-t0)
    
    kwargs['completed']=str(time.time()-t0)
    
    store_parameters(**kwargs)

#-------------------------------------------------------------
if __name__ == "__main__":

    run_model(file_name=project.home+'detail.csv', aa_scenario_name=project.scenario_name,
              ab_time=project.time, res_factor= project.res_factor, tide=project.tide, user=project.user,
              alpha = project.alpha, friction=project.friction,
              time_thinning = project.time_thinning,
              dir_comment=project.dir_comment)