source: anuga_work/production/australia_ph2/melbourne/run_model.py @ 6692

"""Run 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_elevation.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_elevation.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 os.path
import time
from time import localtime, strftime, gmtime

# Related major packages
from Scientific.IO.NetCDF import NetCDFFile
import Numeric as num

from anuga.interface import create_domain_from_regions
from anuga.interface import Transmissive_stage_zero_momentum_boundary
from anuga.interface import Dirichlet_boundary
from anuga.interface import Reflective_boundary
from anuga.interface import Field_boundary
from anuga.interface import Time_boundary
from anuga.interface import file_function

from anuga.interface import create_sts_boundary
from anuga.interface import csv2building_polygons
from file_length import file_length

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 urs2sts
from anuga.utilities.polygon import read_polygon, Polygon_function

# Application specific imports
from setup_model import project
import build_urs_boundary as bub
import prepare_timeboundary as TB

#-------------------------------------------------------------------------------
# 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, __file__, 
                os.path.join(os.path.dirname(project.__file__),
                             project.__name__+'.py'))
start_screen_catcher(project.output_run, 0, 1)

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

# Create the STS file
print 'project.mux_data_folder=%s' % project.mux_data_folder
if not os.path.exists(project.event_sts_east + '.sts'):
    bub.build_urs_boundary(project.mux_input_filename, project.event_sts_east, project.urs_order_east)
if not os.path.exists(project.event_sts_west + '.sts'):
    bub.build_urs_boundary(project.mux_input_filename, project.event_sts_west, project.urs_order_west)

# Read in boundary from ordered sts file
event_sts_east = create_sts_boundary(project.event_sts_east)
event_sts_west = create_sts_boundary(project.event_sts_west)

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

##print 'land N', landward_boundary_N, ' with length ', len(landward_boundary_N)
##print 'land S', landward_boundary_S, ' with length ', len(landward_boundary_S)

# Combine sts polyline with landward points
bounding_polygon_sts = landward_boundary_N + event_sts_east + landward_boundary_S + event_sts_west

print 'bounding polygon sts',bounding_polygon_sts,'length',len(bounding_polygon_sts)

# 4 sides: north (land), east (ocean), south (land), west (ocean)
num_north_points = len(landward_boundary_N)
num_east_points = len(event_sts_east)
num_south_points = len(landward_boundary_S)
num_west_points = len(event_sts_west)

# Boundary tags refer to project.landward_boundary
# 4 points equals 5 segments start at N
boundary_tags={'back': range(num_north_points-1) + range(num_north_points+num_east_points,num_north_points+num_east_points+num_south_points-1),
               'side': [num_north_points-1,
                        num_north_points+num_east_points-1,
                        num_north_points+num_east_points+num_south_points-1,
                        num_north_points+num_east_points+num_south_points+num_west_points-1],
               'oceanE': range(num_north_points,num_north_points+num_east_points-1),
               'oceanW': range(num_north_points+num_east_points+num_south_points,num_north_points+num_east_points+num_south_points+num_west_points-1)}

##print 'tags',boundary_tags

# Build mesh and domain
domain = create_domain_from_regions(bounding_polygon_sts,
                                    boundary_tags=boundary_tags,
                                    maximum_triangle_area=project.bounding_maxarea,
                                    interior_regions=project.interior_regions,
                                    mesh_filename=project.meshes,
                                    use_cache=True,
                                    verbose=True)
print domain.statistics()

domain.set_name(project.scenario_name)
domain.set_datadir(project.output_run) 
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
if project.land_initial_conditions:
    IC = Polygon_function(project.land_initial_conditions,
                          default=project.tide,
                          geo_reference=domain.geo_reference)
else:
    IC = 0
domain.set_quantity('stage', IC, use_cache=True, verbose=True)
domain.set_quantity('friction', project.friction) 
domain.set_quantity('elevation', 
                    filename=project.combined_elevation+'.pts',
                    use_cache=True,
                    verbose=True,
                    alpha=project.alpha)

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

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

# Prepare time boundary
TB.prepare_timeboundary(project.boundary_csv)
f = file_function(project.boundary_csv[:-4] + '.tms')

Br = Reflective_boundary(domain)
Bt = Transmissive_stage_zero_momentum_boundary(domain)
Bd = Dirichlet_boundary([project.tide, 0, 0])

if project.wave == 'Bf':
    Bfe = Field_boundary(project.event_sts_east+'.sts',
                        domain, mean_stage=project.tide,
                        time_thinning=1,
                        default_boundary=Dirichlet_boundary([0, 0, 0]),
                        boundary_polygon=bounding_polygon_sts,                    
                        use_cache=True,
                        verbose=True)
    Bfw = Field_boundary(project.event_sts_west+'.sts',
                        domain, mean_stage=project.tide,
                        time_thinning=1,
                        default_boundary=Dirichlet_boundary([0, 0, 0]),
                        boundary_polygon=bounding_polygon_sts,                    
                        use_cache=True,
                        verbose=True)
    domain.set_boundary({'back': Br,
                     'side': Bt,
                     'oceanE': Bfe,
                     'oceanW': Bfw}) 
    
elif project.wave == 'Tb':
    Tb = Time_boundary(domain,f,default_boundary=Dirichlet_boundary([0, 0, 0]) )

    if project.event_side == 'east' :
        domain.set_boundary({'back': Br,
                             'side': Bt,
                             'oceanW': Bt,
                             'oceanE': Tb})
    elif project.event_side == 'west' :
        domain.set_boundary({'back': Br,
                             'side': Bt,
                             'oceanW': Tb,
                             'oceanE': Bt})
    else:
        print 'need to specify which side event is on'

else:
    print 'No wave specified in project script (Bf or Tb)'
    

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