source: anuga_work/development/attenuation_near_shore/run_beach.py @ 5101

"""

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.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_trial_run=False):
    basename = 'zz' + str(friction)
    if is_trial_run is True:
        outputdir_name += '_test'
        yieldstep = 1
        finaltime = 100.
        maximum_triangle_area=3000
        thinner=True
    else:
        yieldstep = 1.
        finaltime = 1600
        finaltime = 100
        maximum_triangle_area = 100
        thinner=True
        
        
    pro_instance = project.Project(['results'],
                                   outputdir_name=outputdir_name,
                                   home='.')
    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
    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 )
    #-------------------------------------------------------------------------
    # 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(pro_instance.outputdir)
    domain.set_default_order(2) # Use second order spatial scheme
    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])
    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'

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

  
    return pro_instance

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