source: trunk/anuga_work/development/joaqium_luis/run.py @ 7924

"""Script for running a tsunami inundation scenario for Boca do Rio.
Adopted from cairns files

"""

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

# Standard modules
import os
import time
import sys

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water import Transmissive_boundary
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import Time_boundary
from anuga.shallow_water import File_boundary
from anuga.shallow_water import Field_boundary

from anuga.pmesh.mesh_interface import create_mesh_from_regions
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf
from anuga.shallow_water.data_manager import dem2pts
from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
from anuga.abstract_2d_finite_volumes.util import file_function

from anuga.caching import cache
from anuga.utilities.polygon import read_polygon, plot_polygons, \
                                    polygon_area, is_inside_polygon

#------------------------------------------------------------------------------
# Define scenario as either ...
#------------------------------------------------------------------------------
scenario = 'br'

if os.access(scenario, os.F_OK) == 0:
    os.mkdir(scenario)
basename = scenario + 'source'

#------------------
# Initial condition
#------------------
tide = 0.0


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

# Filenames
#dem_name = 'sub_region_bat' 
#meshname = 'sub_region.msh'
dem_name = 'bocarra' 
meshname = 'bocarra.msh'


#---------
# Polygons
#---------

# bounding polygon for study area
bounding_polygon = read_polygon('out_bocarra.csv')

print 'Area of bounding polygon in km?', polygon_area(bounding_polygon)/1000000.0

# interior polygons
#poly_midle = read_polygon('midle_rect.csv')
poly_shallow = read_polygon('shallow.csv')


#------------
# Resolutions
#------------
remainder_res = 2000
midle_res = 3500
shallow_res = 100

interior_regions = [[poly_shallow, shallow_res]]

def setup_domain(tide,
                 bounding_polygon,
                 remainder_res,                 
                 interior_regions,
                 mesh_filename,
                 points_filename):
    """Set up domain except boundary conditions
    This function is defined in order to cache it.
    """





    create_mesh_from_regions(bounding_polygon,
                             boundary_tags={'west': [0],
                                            'north': [1],
                                            'east': [2],
                                            'south': [3]},
                             maximum_triangle_area=remainder_res,
                             filename=mesh_filename,
                             interior_regions=interior_regions,
                             use_cache=False,
                             verbose=True)

    #--------------------------------------------------------------------------
    # Setup computational domain
    #--------------------------------------------------------------------------

    domain = Domain(mesh_filename,
                    use_cache=False,
                    verbose=True)

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

    domain.set_name(basename) 
    domain.set_datadir(scenario)
    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    domain.set_minimum_storable_height(0.05)

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

    domain.set_quantity('stage', tide)
    domain.set_quantity('friction', 0.025) 
    domain.set_quantity('elevation',
                        filename=points_filename,
                        use_cache=False,
                        verbose=True)

    return domain
    

#--------------------------------    
# Call (and cache) setup function
#--------------------------------    

# Run set_up domain without caching
#domain = setup_domain(tide, bounding_polygon, poly_shallow, 
#                      remainder_res, midle_res, shallow_res, 
#                      meshname, 
#                      dem_name + '.pts')

# Run set_up domain with caching
domain = cache(setup_domain,
               (tide,
                bounding_polygon,
                remainder_res, 
                interior_regions, 
                meshname, 
                dem_name + '.pts'),
               #dependencies = 
               verbose=True)            



#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

print 'Available boundary tags', domain.get_boundary_tags()
Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])
Bf = Field_boundary('swan_tsu_stageLand.sww', domain, time_thinning=1, mean_stage=tide,
                    use_cache=True, verbose=True)
   
# Boundary condition for sww feed at the east boundary
domain.set_boundary({'west': Bf,'south':Bf,'east': Br,'north': Br})

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

import time
t0 = time.time()

from Numeric import allclose
from anuga.abstract_2d_finite_volumes.quantity import Quantity

# Save every 1 sec leading up to wave approaching land
for t in domain.evolve(yieldstep = 10, finaltime = 90):
    domain.write_time()
    domain.write_boundary_statistics(tags = ['west','south'])

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