source: anuga_work/production/busselton/standardised_version/run_model.py @ 6368

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



##
# @brief Build boundary STS files from one or more MUX files.
# @param mux_dir Directory containing the MUX files.
# @param event_file Path to meta-file containing mux files+weight data.
# @param output_dir Directory to write STS data to.
# @note 'event_file' is produced by EventSelection.
def build_urs_boundary(mux_dir, event_file, output_dir):
    '''Build a boundary STS file from a set of MUX files.'''

    print 'build_urs_boundary: mux_dir=%s' % mux_dir
    print 'build_urs_boundary: event_file=%s' % event_file
    print 'build_urs_boundary: output_dir=%s' % output_dir

    # if we are using an EventSelection multi-mux file
    if project.multi_mux:    
        # get the mux+weight data from the event file
        mux_event_file = os.path.join(mux_dir, event_file)
        try:
            fd = open(mux_event_file, 'r')
            mux_data = fd.readlines()
            fd.close()
        except IOError, e:
            msg = 'File %s cannot be read: %s' % (mux_event_file, str(e))
            raise Exception, msg
        except:
            raise

        # first line of file is # filenames+weight in rest of file
        num_lines = int(mux_data[0].strip())
        mux_data = mux_data[1:]
        print 'number of sources %d' % num_lines

        # quick sanity check on input mux meta-file
        if num_lines != len(mux_data):
            msg = ('Bad file %s: %d data lines, but line 1 count is %d'
                   % (event_file, len(mux_data), num_lines))
            raise Exception, msg

        # Create filename and weights lists.
        # Must chop GRD filename just after '*.grd'.
        mux_filenames = []
        for line in mux_data:
            muxname = line.strip().split()[0]
            split_index = muxname.index('.grd')
            muxname = muxname[:split_index+len('.grd')]
            muxname = os.path.join(mux_dir, muxname)
            mux_filenames.append(muxname)
        
        mux_weights = [float(line.strip().split()[1]) for line in mux_data]

        # Call legacy function to create STS file.
        # This should be replaced in the future.
        print 'reading', project.urs_order
        print 'creating STS file'
        print 'mux_filenames=%s' % str(mux_filenames)
        print 'basename_out=%s' % str(output_dir)
        print 'ordering_filename=%s' % str(project.urs_order)
        print 'weights=%s' % str(mux_weights)
        print 'mean_stage=%s' % str(project.tide)
        urs2sts(mux_filenames,
                basename_out=output_dir,
                ordering_filename=project.urs_order,
                weights=mux_weights,
                mean_stage=project.tide,
                verbose=True)
    else:                           # a single mux stem file
        urs_filenames = [os.path.join(mux_dir, event_file)]

        weight_factor = 1.0
        weights = weight_factor*num.ones(len(urs_filenames), num.float)
            
        order_filename = os.path.join(project.order_filename_dir)

        print 'reading', order_filename
        # Create ordered sts file
        print 'creating sts file'
        urs2sts(urs_filenames,
                basename_out=os.path.join(project.boundaries_dir,
                                          project.scenario_name),
                ordering_filename=order_filename,
                weights=weights,
                mean_stage=project.tide,
                verbose=True)

    # report on stuff so far
    quantities, elevation, time = get_sts_gauge_data(project.event_folder,
                                                     verbose=False)
    print len(elevation), len(quantities['stage'][0,:])

#-------------------------------------------------------------------------------
# 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
bub.build_urs_boundary(project.mux_input_filename,
                       os.path.join(project.event_folder,
                                    project.scenario_name))

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

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

# Combine sts polyline with landward points
bounding_polygon_sts = event_sts + landward_boundary

# Number of boundary segments
num_ocean_segments = len(event_sts) - 1
# Number of landward_boundary points
num_land_points = file_length(project.landward_boundary)

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

# 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
IC = Polygon_function(project.land_initial_conditions,
                      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_elevation+'.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)
Bt = Transmissive_stage_zero_momentum_boundary(domain)
Bd = Dirichlet_boundary([kwargs['tide'], 0, 0])
Bf = Field_boundary(project.event_sts+'.sts',
                    domain, mean_stage=project.tide,
                    time_thinning=1,
                    default_boundary=Bd,
                    boundary_polygon=bounding_polygon_sts,                    
                    use_cache=True,
                    verbose=True)

domain.set_boundary({'back': Br,
                     'side': Bt,
                     '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)