source: anuga_work/development/near_shore_PMD/run_beach.py @ 5191

"""

Script for running a breaking wave simulation of Jon Hinwoods wave tank.
Note: this is based on the frinction_ua_flume_2006 structure.


Duncan Gray, GA - 2007 



"""


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

# Standard modules
import time
from time import localtime, strftime
import sys
from shutil import copy
from os import path, sep
from os.path import dirname  #, basename
from math import sin, pi

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, \
                            Dirichlet_boundary,  Time_boundary, \
                            File_boundary, \
                            Transmissive_Momentum_Set_Stage_boundary
from anuga.fit_interpolate.interpolate import interpolate_sww2csv
from anuga.abstract_2d_finite_volumes.util import copy_code_files, \
     file_function
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
 

from anuga.shallow_water.data_manager import start_screen_catcher

from anuga.abstract_2d_finite_volumes.generic_boundary_conditions\
     import File_boundary_time

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

def elevation_function(x,y):
    from Numeric import zeros, size, Float

    z = zeros(size(x), Float)
    for i in range(len(x)):
        if x[i] <= 0:
            # swash
            z[i] = -0.05*x[i] 
            #z[i] = -30
        elif x[i] <= 1000:
            # surf
            z[i] = -0.15*(x[i]**(0.5))
            #z[i] = -30
        else:
            # Nominal shoreface
            z[i] = -0.01888*(x[i]**(0.8))
            #z[i] = -30
            #>>> -0.15*(1000**0.5)
            #-4.7434164902525691

            #>>> 0.01888*(1000**0.8)
            #4.7424415826900885

            #>>> -0.01888*10000**0.8
            #-29.922783473665834
    return z 

def main(friction=0.01, outputdir_name=None,
         is_structured=False,
         is_trial_run=False):

    # Structured mesh variables
    dx = 5.
    dy = dx
    L = 2080.
    W = dx
    outputdir_name += '_dx' + str(dx)
    
    basename = 'dx_' + str(dx)
    if is_trial_run is True:
        outputdir_name += '_test'
        yieldstep = 1
        finaltime = 100.
        maximum_triangle_area=3000
        thinner=True
    else:
        yieldstep = 1.
        finaltime = 1600
        maximum_triangle_area = 100
        thinner=True
        
        
    pro_instance = project.Project(['data','near_shore_PMD'],
                                   outputdir_name=outputdir_name)
    print "The output dir is", pro_instance.outputdir
    copy_code_files(pro_instance.outputdir,__file__,
                    dirname(project.__file__) \
                    + sep + project.__name__+'.py')
    copy (pro_instance.codedir + 'run_beach.py',
          pro_instance.outputdir + 'run_beach.py')
    copy (pro_instance.codedir + 'create_mesh.py',
          pro_instance.outputdir + 'create_mesh.py')
    
    #mesh_filename = pro_instance.meshdir + basename + '.tsh'
    mesh_filename = pro_instance.meshdir + basename + '.msh'

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

    # creates copy of code in output dir
    if is_trial_run is False:
        start_screen_catcher(pro_instance.outputdir) #, rank, pypar.size())

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

    # this creates the mesh
    #gate_position = 12.0
    if is_structured is True:

        # Original struc mesh values
        # dx = 10.
        # dy = dx
        # L = 2080.
        # W = dx
        
        # structured mesh
        points, vertices, boundary = rectangular_cross(
            int(L/dx), int(W/dy), L, W, (-80.0, -W*0.5))

        domain = Domain(points, vertices, boundary) 
    else:
        create_mesh.generate(mesh_filename,
                             maximum_triangle_area=maximum_triangle_area,
                             thinner=thinner)

        head,tail = path.split(mesh_filename)
        copy (mesh_filename,
              pro_instance.outputdir + tail )
        domain = Domain(mesh_filename, use_cache = False, verbose = True)
        domain = Domain(points, vertices, boundary) 
    #-------------------------------------------------------------------------
    # Setup computational domain
    #-------------------------------------------------------------------------
   

    print 'Number of triangles = ', len(domain)
    print 'The extent is ', domain.get_extent()
    print domain.statistics()

    
    domain.set_name(basename)
    domain.set_datadir(pro_instance.outputdir)
    domain.set_default_order(2) # Use second order spatial scheme
    domain.set_timestepping_method('rk2')
    domain.set_default_order(2)
    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
    domain.beta_h      = 0.0

    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
    #-------------------------------------------------------------------------

    domain.set_quantity('stage', 0.0)
    domain.set_quantity('friction', friction) 
    domain.set_quantity('elevation', elevation_function)

    
    print 'Available boundary tags', domain.get_boundary_tags()

    # Create boundary function from timeseries provided in file
    #function = file_function(project.boundary_file, domain, verbose=True)
    #Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
    Br = Reflective_boundary(domain)
    Bd = Dirichlet_boundary([10.,0,0])
    Bwp = Transmissive_Momentum_Set_Stage_boundary(domain,
                   lambda t: [(1.2*sin(2*t*pi/10)), 0.0 ,0.0])
    Bwp_velocity = Time_boundary(domain,
                   lambda t: [(1.2*sin(2*t*pi/10)),
                              0.24*sin(2*t*pi/10),0.0])
    Bws = Transmissive_Momentum_Set_Stage_boundary(domain,
                   lambda t: [1.2*sin(2*t*pi/10)+1.0*sin(2*t*pi/9.5),
                              0.0, 0.0])
    
    if is_structured is True:
        domain.set_boundary( {'top': Br,
                              'bottom': Br ,
                              'left' : Br ,
                              'right' : Bwp } )
    else:
        domain.set_boundary( {'wall': Br, 'wave': Bwp} )


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

    points = [[-2,0.0],  #-1.8m from SWL
              [0,0.0],  #0.5m from SWL
              [2,0.0],  #2m from SWL
              ]


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

    if True:
        interpolate_sww2csv(pro_instance.outputdir + basename +".sww",
                            points,
                            pro_instance.outputdir + "depth_dx_"+str(dx)+".csv",
                            pro_instance.outputdir + "velocity_x_dx_"+str(dx)+".csv",
                            pro_instance.outputdir + "velocity_y_dx_"+str(dx)+".csv")
  
    return pro_instance

#-------------------------------------------------------------
if __name__ == "__main__":
    main(
        is_trial_run=False,
        friction=0.0,
        is_structured=True,
        outputdir_name='structured_old_limiters')