source: trunk/anuga_work/development/uqwave/run_UQwave.py

"""Script for running a simulation of UQ's wave flume.


Ole Nielsen and Duncan Gray, GA - 2006 
Modified by Matt Barnes 2012
"""


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

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



import anuga
from anuga.abstract_2d_finite_volumes.region import Set_region  
from anuga.fit_interpolate.interpolate import interpolate_sww2csv, \
    file_function
from anuga.utilities.file_utils import copy_code_files
from anuga.file.csv_file import load_csv_as_dict
from numerical_tools import  err
 
from anuga.file_conversion.file_conversion import timefile2netcdf


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


def main():

    D = 0.2 #initial depth
    frictions = [0.009]
    
    for friction in frictions:
        scenario(D, friction)
                        
def scenario(D, friction):

    import project
    #-------------------------------------------------------------------------
    # Setup archiving of simulations
    #-------------------------------------------------------------------------

    id = 'D_'+str(D)+'f_'+str(friction)
    copy (project.codedirname, project.outputtimedir + 'project.py')
    run_name = 'run_dam.py'
    run_name_out = 'run_dam'+id+'.py'
    
    copy (project.codedirname, project.outputtimedir + 'project.py')
    copy (project.codedir + 'run_UQwave.py', project.outputtimedir + 'run_UQwave.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
    #-------------------------------------------------------------------------
    flume_length = 7.24
    #flume_length = 3.24


    create_mesh.generate(project.mesh_filename,
                         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 = anuga.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':2,
                                        'xmomentum': 2,
                                        'ymomentum': 2,
                                        'elevation': 1})
    domain.set_minimum_storable_height(0.003)
    domain.set_store_vertices_uniquely(True)  # for writting to sww

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

    slope = 0.084  #beach slope
    
    def elevation_tilt(x, y):
        return x*slope
        
    domain.set_quantity('stage', D)    #water level
    domain.set_quantity('friction', friction) #bed friction
    domain.set_quantity('elevation', 0)

    domain.set_region('beach', 'elevation', elevation_tilt)
                                      
    print 'Available boundary tags', domain.get_boundary_tags()
   
    Br = anuga.Reflective_boundary(domain)
    Bd = anuga.Dirichlet_boundary([0,0,0])  # to drain the water out.

    #Bore file input

    bore_file = 'bore'
    timefile2netcdf(bore_file + '.txt',
                    quantity_names = ['stage','xmomentum', 'ymomentum'],
                    time_as_seconds=True)
    F = anuga.File_boundary(bore_file + '.tms', domain)
    
    domain.set_boundary( {'wall': Br, 'back':F} )

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

    for t in domain.evolve(yieldstep = 0.02, finaltime = 10): #enter timestep and final time
        domain.write_time()
    
        print 'That took %.2f seconds' %(time.time()-t0)
        print 'finished'

    points = []
    for i in range(41):
        points.append([float(i)/20 + -2.05, 0.425])

    for j in range(120):
        points.append([float(j)/20, 0.425])
              
    #-------------------------------------------------------------------------
    # Calculate gauge info
    #-------------------------------------------------------------------------
    if True:
        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()