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

"""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, sep
from os.path import dirname  #, basename

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, \
                            Dirichlet_boundary, Time_boundary, File_boundary
from anuga.abstract_2d_finite_volumes.region import Set_region 
from anuga.fit_interpolate.interpolate import interpolate_sww2csv
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, \
     copy_code_files


# Scenario specific imports
import project                 # Definition of file names and polygons
import create_mesh


def main(friction):

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

    # creates copy of code in output dir

    print "The output dir is", project.outputtimedir
    copy_code_files(project.outputtimedir,__file__,
                    dirname(project.__file__)+sep+ project.__name__+'.py' )
    copy (project.codedir + 'create_mesh.py',
          project.outputtimedir + 'create_mesh.py')
    myid = 0
    numprocs = 1
    #start_screen_catcher(project.outputtimedir, myid, numprocs)

    print 'USER:    ', project.user
    #-------------------------------------------------------------------------
    # Create the triangular mesh
    #-------------------------------------------------------------------------

    gate_position = 0.85
    create_mesh.generate(project.mesh_filename,
                         gate_position,
                         #is_coarse=True) # this creates the mesh
                         is_coarse=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_storable_height(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', friction) 
    domain.set_quantity('elevation',elevation_tilt)

    print 'Available boundary tags', domain.get_boundary_tags()
    domain.set_region('dam','stage',0.20,
                                 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.01, finaltime = 30):
        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]
              ]

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