source: anuga_work/development/gong_2008/run_gong_slide.py @ 6868

"""Script for running a tsunami inundation scenario for Wollongong, NSW, Australia.

Source data such as elevation and boundary data is assumed to be available in
directories specified by project_slide.py
The output sww file is stored in project_slide.outputtimedir 

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and a tsunami wave generated by s submarine mass failure.

Ole Nielsen and Duncan Gray, GA - 2005 and Nick Bartzis, GA - 2006
"""

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

# Standard modules
import os
import time
from shutil import copy
from os.path import dirname, basename
from os import mkdir, access, F_OK, sep
import sys

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
from anuga.shallow_water.data_manager import convert_dem_from_ascii2netcdf, dem2pts
from anuga.geospatial_data.geospatial_data import *
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files

# Application specific imports
import project_slide              # Definition of file names and polygons

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


#-------------------------------------------------------------------------------
# Preparation of topographic data
#
# Convert ASC 2 DEM 2 PTS using source data and store result in source data
#-------------------------------------------------------------------------------
"""
# filenames
on_offshore10_dem_name = project_slide.on_offshore10_dem_name
nsw_dem_name = project_slide.nsw_dem_name


# creates DEM from asc data
convert_dem_from_ascii2netcdf(on_offshore10_dem_name, use_cache=True, verbose=True)
convert_dem_from_ascii2netcdf(nsw_dem_name, use_cache=True, verbose=True)

#creates pts file for onshore DEM
dem2pts(on_offshore10_dem_name, use_cache=True, verbose=True)
dem2pts(nsw_dem_name,
        easting_min=project_slide.eastingmin_nsw,
        easting_max=project_slide.eastingmax_nsw,
        northing_min=project_slide.northingmin_nsw,
        northing_max= project_slide.northingmax_nsw,
        use_cache=True, verbose=True)

print 'create offshore'
G11 = Geospatial_data(file_name = project_slide.offshore_dem_name1 + '.txt')
G12 = Geospatial_data(file_name = project_slide.offshore_dem_name4 + '.txt')+\
      Geospatial_data(file_name = project_slide.offshore_dem_name5 + '.txt')+\
      Geospatial_data(file_name = project_slide.offshore_dem_name6 + '.txt')+\
      Geospatial_data(file_name = project_slide.offshore_dem_name7 + '.txt')+\
      Geospatial_data(file_name = project_slide.offshore_dem_name8 + '.txt')+\
      Geospatial_data(file_name = project_slide.offshore_dem_name9 + '.txt')
print 'create onshore'
G2 = Geospatial_data(file_name = project_slide.on_offshore10_dem_name + '.pts')
G4 = Geospatial_data(file_name = project_slide.nsw_dem_name + '.pts')
print 'add'
G = G11.clip(Geospatial_data(project_slide.poly_surveyclip)) +\
    G12.clip(Geospatial_data(project_slide.polyAll)) +\
    G2.clip(Geospatial_data(project_slide.poly_10mclip)) +\
    (G4.clip(Geospatial_data(project_slide.polyAll))).clip_outside(Geospatial_data(project_slide.poly_surveyclip)).clip_outside(Geospatial_data(project_slide.poly_10mclip))
print 'export points'
G.export_points_file(project_slide.combined_dem_name + '.pts')
#G.export_points_file(project_slide.combined_dem_name + '.xya')
"""
#----------------------------------------------------------------------------
# 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
#-------------------------------------------------------------------------------

from anuga.pmesh.mesh_interface import create_mesh_from_regions
meshname = project_slide.meshname+'.msh'
remainder_res = 100000*100
local_res = 25000*1000
gong_res = 500*1000
interior_regions = [[project_slide.poly_local, local_res],
                    [project_slide.poly_gong, gong_res],
                    [project_slide.poly_southgong, gong_res]]

from caching import cache
_ = cache(create_mesh_from_regions,
          project_slide.polyAll,
           {'boundary_tags': {'e0': [0], 'e1': [1], 'e2': [2],
                              'e3': [3], 'e4':[4]},
           'maximum_triangle_area': remainder_res,
           'filename': meshname,
           'interior_regions': interior_regions},
          verbose = True, evaluate=False)
print 'created mesh'

#-------------------------------------------------------------------------------                                 
# 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(project_slide.basename)
domain.set_datadir(project_slide.outputtimedir)
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)

#-------------------------------------------------------------------------------                                 
# Set elevation to mesh
#-------------------------------------------------------------------------------

tide = 0.0
domain.set_quantity('elevation', 
                    filename = project_slide.combined_dem_name + '.pts',
                    use_cache = True,
                    verbose = True,
                    alpha = 0.1
                    )

#-------------------------------------------------------------------------------                                 
# Setup boundary conditions
#-------------------------------------------------------------------------------
print 'Available boundary tags', domain.get_boundary_tags()

Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([tide,0,0])

domain.set_boundary( {'e0': Bd,  'e1': Bd, 'e2': Bd, 'e3': Bd, 'e4': Bd} )


#-------------------------------------------------------------------------------
# Set up scenario (tsunami_source is a callable object used with set_quantity)
#-------------------------------------------------------------------------------
from smf import slide_tsunami

# effect on a3D and wavelength
l0 = project_slide.birubi_length
w0 = project_slide.birubi_width
d0 = project_slide.birubi_depth
s0 = project_slide.birubi_slope
t0 = project_slide.birubi_thickness
a0 = project_slide.birubi_alpha
x0 = project_slide.slide_origin_birubi_a[0]
y0 = project_slide.slide_origin_birubi_a[1]
length = l0
width = w0
depth = d0
slope = s0
thickness = t0
alpha = a0
gamma0 = 1.85
gamma = gamma0
m0 = 1.
massco = m0
drag0 = 1.
dragco = drag0

# no effect on a3D and wavelength but used in Double Gaussian
dx = 0.01
kappa0 = 3.
kappad0 = 0.8
kappa = kappa0
kappad = kappad0

# this doesn't seem to apper anywhere in smf
frictionc0 = 0.01
frictionco = frictionc0

# scaling for Double Gaussian function
scale0 = 100. # Bridgette's fiddle 
scale = scale0 + 10.

for i in range(1):
    scale = scale - 10.   
    mydir = project_slide.outputdir+'testing_' + str(int(scale))
    print 'dir_comment', mydir
    start_screen_catcher(mydir, 0, 1)
    # creates copy of code in output dir 
    copy_code_files(mydir,__file__,dirname(project_slide.__file__)+sep+ project_slide.__name__+'.py' )
    start_screen_catcher(mydir, 0, 1)

    print 'USER:    ', project_slide.user

    tsunami_source = slide_tsunami(length=length,
                                   width=width,
                                   depth=depth,
                                   slope=slope,
                                   thickness=thickness, 
                                   x0=x0, 
                                   y0=y0, 
                                   alpha=alpha,
                                   gamma=gamma,
                                   massco=massco,
                                   dragco=dragco,
                                   frictionco=frictionco,
                                   dx=dx,
                                   kappa=kappa,
                                   kappad=kappad,
                                   scale=scale,
                                   domain=domain,
                                   verbose=True)

    print 'hello', scale, tsunami_source.wavelength, tsunami_source.a3D
    
    #-------------------------------------------------------------------------------                                 
    # Setup initial conditions
    #-------------------------------------------------------------------------------

    domain.set_quantity('stage', tsunami_source)
    domain.set_quantity('friction', 0.0) 


    #-------------------------------------------------------------------------------                                 
    # Evolve system through time
    #-------------------------------------------------------------------------------
    import time
    t0 = time.time()
    from Numeric import allclose
    from anuga.abstract_2d_finite_volumes.quantity import Quantity

    for t in domain.evolve(yieldstep = 30, finaltime = 5000): 
        domain.write_time()
        domain.write_boundary_statistics(tags = 'e14')
        
        
    print 'That took %.2f seconds' %(time.time()-t0)

    print 'finished'