source: anuga_work/development/joaqium_luis/run.py @ 4925

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


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

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

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

# Create DEM from asc data
#convert_dem_from_ascii2netcdf(dem_name, use_cache=True, verbose=True)

# Create pts file for onshore DEM
#dem2pts(dem_name, use_cache=True, verbose=True)

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

###############################
# Domain definitions
###############################

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

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

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

remainder_res = 2000
midle_res = 3500
shallow_res = 100
interior_regions = [[bounding_polygon, remainder_res], [poly_shallow, shallow_res]]
#interior_regions = [[bounding_polygon, remainder_res], [poly_midle, midle_res], [poly_shallow, shallow_res]]

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

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

domain = Domain(meshname, use_cache=True, 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_quantities_to_be_stored('stage')
domain.set_minimum_storable_height(0.05)

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

tide = 0.0
domain.set_quantity('stage', tide)
domain.set_quantity('friction', 0.025) 
domain.set_quantity('elevation', filename=dem_name + '.pts',
                    use_cache=True, verbose=True, alpha=0.1)

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

# Create boundary function from timeseries provided in file
#function1 = file_function('mareg.tms', domain, verbose=True)

# Create and assign boundary objects
#Bts1 = Transmissive_Momentum_Set_Stage_boundary(domain, function1) 

print 'Available boundary tags', domain.get_boundary_tags()
#Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
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)