source: anuga_validation/Hinwood_2008/run_dam.py @ 5689

"""

Script for running a breaking wave simulation of Jon Hinwoods wave tank.

Duncan Gray, GA - 2007 
"""


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

# Standard modules
import time


# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, Time_boundary
from anuga.fit_interpolate.interpolate import interpolate_sww2csv
from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.utilities.interp import interp

# Scenario specific imports
import create_mesh
from prepare_time_boundary import csv2tms


def main(metadata_dic,
         maximum_triangle_area,
         yieldstep,
         boundary_path=None,
         friction=0.,  # planed wood 0.012 http://www.lmnoeng.com/manningn.htm
         outputdir_name=None,
         isTest=False,
         width=1.0,
         use_limits=True,
         end_tag = ''):

    id = metadata_dic['scenario_id']
    
    if isTest is True:
        yieldstep = 1.0
        finaltime = 15.
        maximum_triangle_area=0.1
        outputdir_name += '_test'         
    else:
        finaltime = None
        
    outputdir_name = paras2outputdir_tag(mta=maximum_triangle_area,
                        yieldstep=yieldstep,
                        width=width,
                        use_limits=use_limits,
                        friction=friction,
                        end_tag=end_tag,
                        outputdir_name=outputdir_name)
    basename = outputdir_name    
    metadata_dic = set_z_origin_to_water_depth(metadata_dic)    

    if finaltime is None:
        finaltime = metadata_dic['wave_times'][1] + 0.1 
        
    boundary_file_path = id + '_boundary.tsm'
        
    # Convert the boundary file, .csv to .tsm
    csv2tms(boundary_file_path, metadata_dic['xleft'][1])
   
    mesh_filename =  basename + '.msh'

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

    create_mesh.generate(mesh_filename, metadata_dic, width=width,
                         maximum_triangle_area=maximum_triangle_area)
    
    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
    domain = Domain(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(basename)
    domain.set_datadir('.')
    domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
    domain.set_minimum_storable_height(0.0001)

    if use_limits is True:
        domain.set_default_order(2) # Use second order spatial scheme
        domain.set_timestepping_method('rk2')
        domain.use_edge_limiter = True
        domain.tight_slope_limiters = True
        
        domain.beta_w      = 0.6
        domain.beta_uh     = 0.6
        domain.beta_vh     = 0.6
    

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

    domain.set_quantity('stage', 0.) #the origin is the still water level
    domain.set_quantity('friction', friction)
    elevation_function = Elevation_function(metadata_dic)
    domain.set_quantity('elevation', elevation_function)

    function = file_function(boundary_file_path, domain, verbose=True)
        
    Br = Reflective_boundary(domain)
    Bts = Time_boundary(domain, function)
    domain.set_boundary( {'wall': Br, 'wave': Bts} )

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

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


    #-------------------------------------------------------------------------
    # Calculate gauge info
    #-------------------------------------------------------------------------

    if isTest is not True:
        points = []
        for gauge_x in metadata_dic['gauge_x']:
            points.append([gauge_x, 0.0])

        id =  ".csv"
        interpolate_sww2csv(basename +".sww",
                            points,
                            basename + "_depth" + id,
                            basename + "_velocity_x" + id,
                            basename + "_velocity_y" + id,
                            basename + "_stage" + id)
  

def paras2outputdir_tag(mta,
                        yieldstep,
                        width,
                        use_limits,
                        friction,
                        end_tag,
                        outputdir_name=None):
    """
    
    """
    if outputdir_name is None:
        outputdir_name = ''
        
    if use_limits is True:
        outputdir_name += '_lmts'
    else:
        outputdir_name += '_nolmts'
    outputdir_name += '_wdth_' + str(width)
    outputdir_name += '_z_' + str(friction)
    outputdir_name += '_ys_' + str(yieldstep)
    outputdir_name += '_mta_' + str(mta)
    outputdir_name += end_tag

    return outputdir_name


class Elevation_function:
    """
    """
    def __init__(self, slope):
        self.xslope_position = [slope['xleft'][0],slope['xtoe'][0],
                  slope['xbeach'][0],slope['xright'][0]]
        self.yslope_height = [slope['xleft'][1],slope['xtoe'][1],
                  slope['xbeach'][1],slope['xright'][1]]
        
    def __call__(self, x,y):       
        z = interp(self.yslope_height, self.xslope_position, x)
        return z 
    
def set_z_origin_to_water_depth(metadata_dic):
    """
    """
    offset = metadata_dic['offshore_water_depth']
    keys = ['xleft', 'xtoe', 'xbeach', 'xright']
    for x in keys:
            metadata_dic[x][1] -= offset
    return metadata_dic

#-------------------------------------------------------------
if __name__ == "__main__":
    
    from scenarios import scenarios

    #scenarios = [scenarios[0]]
    
    width = 0.1
    maximum_triangle_area=0.01
    yieldstep = 0.01
    #yieldstep = 0.5
    friction=0.0
    isTest=True
    isTest=False
    
    for run_data in scenarios:
        main(run_data,
             maximum_triangle_area=maximum_triangle_area,
             yieldstep=yieldstep,
             width=width,
             isTest=isTest,
             outputdir_name=run_data['scenario_id'],
             use_limits=False,
             friction=friction,
             end_tag='')