source: development/dam_2006/run_dam.py @ 3456

"""Script for running a tsunami inundation scenario for Onslow, WA, Australia.

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

The scenario is defined by a triangular mesh created from project.polygon,
the elevation data and a simulated submarine landslide.

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


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

# Standard modules
import time
import sys
from shutil import copy
from os import mkdir, access, F_OK, path

# Related major packages
from pyvolution.shallow_water import Domain, Reflective_boundary, \
                            Dirichlet_boundary, Time_boundary, File_boundary
from pyvolution.util import Screen_Catcher
from pyvolution.region import Set_region 
from fit_interpolate.interpolate import interpolate_sww2csv
import create_mesh

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

def main():
    #-------------------------------------------------------------------------
    # Setup archiving of simulations
    #-------------------------------------------------------------------------

    copy (project.codedirname, project.outputtimedir + 'project.py')
    copy (project.codedir + 'run_dam.py', project.outputtimedir + 'run_dam.py')
    copy (project.codedir + 'create_mesh.py',
          project.outputtimedir + 'create_mesh.py')
    print'output dir', project.outputtimedir

    #FIXME this isn't working
    #normal screen output is stored in 
    screen_output_name = project.outputtimedir + "screen_output.txt"
    screen_error_name = project.outputtimedir + "screen_error.txt"

    #-------------------------------------------------------------------------
    # Create the triangular mesh
    #-------------------------------------------------------------------------

    gate_position = 0.85
    create_mesh.generate(project.mesh_filename,
                         gate_position,
                         is_course=True) # this creates the mesh
                         #is_course=False) # this creates the mesh

    head,tail = path.split(project.mesh_filename)
    copy (project.mesh_filename,
          project.outputtimedir + tail )
    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
    domain = Domain(project.mesh_filename, use_cache = False, verbose = True)
   

    print 'Number of triangles = ', len(domain)
    print 'The extent is ', domain.get_extent()
    print domain.statistics()

    domain.set_name(project.basename)
    domain.set_datadir(project.outputtimedir)
    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])

    #-------------------------------------------------------------------------
    # Setup initial conditions
    #-------------------------------------------------------------------------

    tide = 0.0
    slope = 0.05

    def elevation_tilt(x, y):
        return x*slope
        
    domain.set_quantity('stage', elevation_tilt)
    domain.set_quantity('friction', 0.03) 
    domain.set_quantity('elevation',elevation_tilt)

    print 'Available boundary tags', domain.get_boundary_tags()
    domain.set_region('dam','stage',0.2,
                                 location = 'unique vertices') 
    domain.set_region(Set_region('dam','stage',0.2,
                                 location = 'unique vertices'))
    domain.set_store_vertices_uniquely(True)  # for writting to sww

    Br = Reflective_boundary(domain)
    Bd = Dirichlet_boundary([0,0,0])
    domain.set_boundary( {'wall': Br, 'edge': Bd} )

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

    for t in domain.evolve(yieldstep = 0.1, finaltime = 25):
        domain.write_time()
    
        print 'That took %.2f seconds' %(time.time()-t0)
        print 'finished'

    points = [[gate_position - 0.65,0.2],
              [gate_position - 0.55,0.2],
              [gate_position - 0.45,0.2],
              [gate_position - 0.35,0.2],
              [gate_position - 0.25,0.2],
              [0.55,0.2]
              ]

    #-------------------------------------------------------------------------
    # Calc gauge info
    #-------------------------------------------------------------------------
    interpolate_sww2csv(project.outputtimedir + project.basename +".sww",
                        points,
                        project.depth_filename,
                        project.velocity_x_filename,
                        project.velocity_y_filename)

#-------------------------------------------------------------
if __name__ == "__main__":
    main()