source: anuga_work/production/busselton/standardised_version/run_busselton.py @ 6262

"""Script for running a tsunami inundation scenario for busselton, 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_busselton.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_busselton.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
import os
import time

# Related major packages
from anuga.interface import create_domain_from_regions
from anuga.interface import Dirichlet_boundary
from anuga.interface import Reflective_boundary
from anuga.interface import Field_boundary
from anuga.interface import create_sts_boundary
from anuga.interface import csv2building_polygons

from anuga.shallow_water.data_manager import start_screen_catcher
from anuga.shallow_water.data_manager import copy_code_files
from anuga.utilities.polygon import read_polygon, Polygon_function
    
# Application specific imports
import project  # Definition of file names and polygons


#-----------------------------------------------------------------------
# Copy scripts to time stamped output directory and capture screen
# output to file
#-----------------------------------------------------------------------

copy script must be before screen_catcher
copy_code_files(project.output_run_time_dir, __file__, 
                os.path.dirname(project.__file__)+os.sep+ project.__name__+'.py' )
start_screen_catcher(project.output_run_time_dir, myid, numprocs)


#--------------------------------------------------------------------------
# Create the computational domain 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
#--------------------------------------------------------------------------
print 'Create computational domain'

# Read in boundary from ordered sts file
urs_bounding_polygon = create_sts_boundary(project.urs_boundary_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

# Number of boundary 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,N+4, N+5],
               'side': [N,N+6],
               'ocean': range(N)}

# Build mesh and domain
domain = create_domain_from_regions(bounding_polygon,
                                    boundary_tags=boundary_tags,
                                    maximum_triangle_area=project.res_poly_all,
                                    interior_regions=project.interior_regions,
                                    mesh_filename=project.meshes_dir_name,
                                    use_cache=True,
                                    verbose=True)
print domain.statistics()

domain.set_name(project.scenario_name)
domain.set_datadir(project.output_run_time_dir) 
domain.set_minimum_storable_height(0.01)    # Don't store depth less than 1cm

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

# Set the initial stage in the offcoast region only
IC = Polygon_function([(project.poly_mainland, 0),
                       (project.poly_marina, 0)], 
                      default=project.tide,
                      geo_reference=domain.geo_reference)
domain.set_quantity('stage', IC, use_cache=True, verbose=True)
domain.set_quantity('friction', project.friction) 
domain.set_quantity('elevation', 
                    filename=project.combined_dir_name+'.pts',
                    use_cache=True,
                    verbose=True,
                    alpha = project.alpha)


#-------------------------------------------------------------------------
# Setup boundary conditions 
#-------------------------------------------------------------------------
print 'Set boundary - available tags:', domain.get_boundary_tags()

Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([project.tide,0,0])
Bf = Field_boundary(project.urs_boundary_name+'.sts',
                    domain, mean_stage=project.tide,
                    time_thinning=1,
                    default_boundary=Bd,
                    boundary_polygon=bounding_polygon,                    
                    use_cache=True,
                    verbose=True)

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


#-------------------------------------------------------------------------
# Evolve system through time
#-------------------------------------------------------------------------
t0 = time.time()
for t in domain.evolve(yieldstep=project.yieldstep, 
                       finaltime=project.finaltime,
                       skip_initial_step=False): 
    print domain.timestepping_statistics()
    print domain.boundary_statistics(tags='ocean')

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