source: anuga_validation/UQ_runup_2006/run_dam.py @ 4678

"""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.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
import create_mesh

# Application specific imports


def main():
     
    slopes = [0]
    frictions = [0.005, 0.01, 0.03] #can enter multiple frictions eg [0.00, 0.01]
    inital_depths = [0.2] #can enter multiple initial depths eg [0.1, 0.2]
    gate_positions = [0.75]
    flume_lengths = [3] #can enter multiple flume lengths eg [3.0, 5.0]
    
    for slope in slopes:
        for friction in frictions:
            for inital_depth in inital_depths:
                for flume_length in flume_lengths:
                    for gate_position in gate_positions:
                        scenario(slope, friction, inital_depth, gate_position, flume_length)
            

def scenario(slope, friction, inital_depth, gate_position, flume_length):
    
    import project                 # Definition of file names and polygons

    #-------------------------------------------------------------------------
    # Setup archiving of simulations
    #-------------------------------------------------------------------------

    id = 'd'+ str(inital_depth) + '_s'+str(slope)+'_g'+ str(gate_position)+ '_n'+ str(friction)
    copy (project.codedirname, project.outputtimedir + 'project.py')
    run_name = 'run_dam.py'
    run_name_out = 'run_dam'+id+'.py'

    copy (project.codedir + run_name, project.outputtimedir + run_name_out)
    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
    #-------------------------------------------------------------------------

    create_mesh.generate(project.mesh_filename,
                         gate_position, flume_length,
                         #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 + id)
    domain.set_datadir(project.outputtimedir)
    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    domain.set_minimum_storable_height(0.001)
    domain.set_store_vertices_uniquely(True)  # for writting to sww

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


    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',inital_depth,
                                 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 = 10): #enter timestep and final time
        domain.write_time()
    
        print 'That took %.2f seconds' %(time.time()-t0)
        print 'finished'

    points = [[0.45, 0.2]]
              

    #-------------------------------------------------------------------------
    # Calculate gauge info
    #-------------------------------------------------------------------------
    #interpolate_sww2csv(project.outputtimedir + project.basename \
     #              + id+".sww",
      #                  points,
       #                 project.depth_filename + id + '.csv',
        #                project.velocity_x_filename + id + '.csv',
         #               project.velocity_y_filename + id + '.csv')

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