source: anuga_work/production/patong/run_patong.py @ 6748

"""Script for running a tsunami inundation scenario for Patong Beach, Thailand.

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data is compiled into a pts file through build_patong.py
and a simulated tsunami for the 2004 event is generated through 
an sts file from build_boundary.py.

Input: sts file (build_boundary.py)
       pts file (build_patong.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

"""

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

# Standard modules
from os import sep
import os
from os.path import dirname
import time

# Related major packages
from anuga.interface import create_domain_from_regions
from anuga.interface import Dirichlet_boundary
from anuga.interface import Transmissive_stage_zero_momentum_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 file_length import file_length

from anuga.caching import cache
from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters
from anuga.utilities.polygon import read_polygon
from polygon import Polygon_function
    
# Application specific imports
import project  # Definition of file names and polygons
numprocs = 1
myid = 0


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

#copy script must be before screen_catcher

output_dir = project.output_run_time_dir
print 'output_dir', output_dir

copy_code_files(output_dir,__file__, 
                dirname(project.__file__)+sep+ project.__name__+'.py' )

start_screen_catcher(output_dir, myid, numprocs)


#-----------------------------------------------------------------------
# Domain definitions
#-----------------------------------------------------------------------

# Read in boundary from ordered sts file
urs_boundary_name = os.path.join(project.boundaries_dir,
                                 project.scenario_name)
urs_bounding_polygon = create_sts_boundary(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

# Counting segments
N = len(urs_bounding_polygon)-1
# Number of landward_boundary points
M = file_length(project.landward_dir)


# Boundary tags refer to project.landward_boundary
# 4 points equals 5 segments start at N
boundary_tags={'back': range(N+1,N+M),
               'side': [N,N+M],
               'ocean': range(N)}

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

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

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

# sets the initial stage in the offcoast region only
IC = Polygon_function([(project.poly_mainland, 0)],
                      default=project.tide,
                      geo_reference=domain.geo_reference)
print 'stage'
domain.set_quantity('stage', IC,
                    use_cache=True,
                    verbose=True)
print 'friction'
domain.set_quantity('friction', project.friction)

print 'elevation'
domain.set_quantity('elevation',
                    filename=project.combined_dir_name+'.pts',
                    alpha=project.alpha,                    
                    use_cache=True,
                    verbose=True)

if project.use_buildings:
    # Add buildings from file
    print 'Reading building polygons'    
    building_polygons, building_heights = csv2building_polygons(project.building_polygon_file)
    #clipping_polygons=project.building_area_polygons)

    print 'Creating %d building polygons' % len(building_polygons)
    def create_polygon_function(building_polygons, geo_reference=None):
        L = []
        for i, key in enumerate(building_polygons):
            if i%100==0: print i
            poly = building_polygons[key]
            elev = building_heights[key]
            L.append((poly, elev))
            
            buildings = Polygon_function(L, default=0.0,
                                         geo_reference=geo_reference)
        return buildings

    buildings = cache(create_polygon_function,
                      building_polygons,
                      {'geo_reference': domain.geo_reference},
                      verbose=True)

    print 'Adding buildings'
    domain.add_quantity('elevation',
                        buildings,
                        use_cache=True,
                        verbose=True)



#------------------------------------------------------
# Distribute domain to implement parallelism !!!
#------------------------------------------------------

#  if numprocs > 1:
#      domain=distribute(domain)

#------------------------------------------------------
# Set domain parameters
#------------------------------------------------------
domain.set_name(project.scenario_name)
domain.set_datadir(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

#-------------------------------------------------------------------------
# Setup boundary conditions 
#-------------------------------------------------------------------------
print 'Set boundary conditions'

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

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


#----------------------------------------------------------------------------
# Evolve system through time
#--------------------------------------------------------------------
t0 = time.time()

# Skip over the first 6000 seconds
for t in domain.evolve(yieldstep=2000,
                       finaltime=6000):
    print domain.timestepping_statistics()
    print domain.boundary_statistics(tags='ocean')

# Start detailed model
for t in domain.evolve(yieldstep=project.yieldstep,
                       finaltime=project.finaltime,
                       skip_initial_step=True):
    print domain.timestepping_statistics()
    print domain.boundary_statistics(tags='ocean')
    
print 'Simulation took %.2f seconds' %(time.time()-t0)