source: anuga_work/development/dam_2006/run_dam.py @ 3610

"""Script for running a dam break simulation of UQ's dam break tank.


Ole Nielsen and Duncan Gray, GA - 2006 
"""


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

# Standard modules
import time
import sys
from shutil import copy
from os import path

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, \
                            Dirichlet_boundary, Time_boundary, File_boundary
from anuga.abstract_2d_finite_volumes.util import Screen_Catcher
from anuga.abstract_2d_finite_volumes.region import Set_region 
from anuga.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'])
    domain.set_minimum_sww_depth(0.01)
    domain.set_store_vertices_uniquely(True)  # for writting to sww

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

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

    Br = Reflective_boundary(domain)
    Bd = Dirichlet_boundary([0,0,0])  # to drain the water out.
    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]
              ]

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