"""Script for running the 3rd Long Wave Workshop Benchmark Problem 4.

Bridgette Lewis and Jane Sexton, July 2008
"""

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

# Standard modules
import os
import time
from shutil import copy
from os.path import dirname, basename
from os import mkdir, access, F_OK, sep
import sys
import project

# Related major packages
from anuga.shallow_water import Domain, Reflective_boundary, Dirichlet_boundary 
from anuga.abstract_2d_finite_volumes.util import start_screen_catcher, copy_code_files
from anuga.shallow_water.data_manager import start_screen_catcher, copy_code_files,store_parameters

myid=0
numprocs=1 
print 'output_dir',project.output_run_time_dir
copy_code_files(project.output_run_time_dir,'benchmark_problem4.py', 
             'project.py' )
start_screen_catcher(project.output_run_time_dir, myid, numprocs)

#----------------------------------------------------------------------------
# Create the triangular mesh based on the problem statement
#-------------------------------------------------------------------------------

from anuga.pmesh.mesh_interface import create_mesh_from_regions
remainder_res = 250000.
#local_res = 50000.
#coast_res = 500.
#interior_regions = [[project_slide.poly_syd1, local_res],
#                    [project_slide.poly_coast, coast_res]]
length = 104. #wave tank length
width = 3.7 #wave tank width

bounding_polygon = [[0,0],[length,0],[length,width],[0,width]]
meshname = 'slide.msh'


from caching import cache
_ = cache(create_mesh_from_regions,
          bounding_polygon,
           {'boundary_tags': {'bottom': [0], 'right': [1], 
                              'top': [2], 'left': [3]},
           'maximum_triangle_area': project.remainder_res,
           'filename': meshname},
           #'interior_regions': interior_regions},
          verbose = True, evaluate=False)
print 'created mesh'

#-------------------------------------------------------------------------------                                 
# Setup computational domain
#-------------------------------------------------------------------------------                                 
domain = Domain(meshname, use_cache = True, verbose = True)

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

domain.set_datadir(project.output_run_time_dir)
domain.set_name(project.sww_name)
domain.set_quantities_to_be_stored(['stage', 'xmomentum', 'ymomentum'])
domain.set_minimum_storable_height(0.01)

#-------------------------------------------------------------------------------                                 
# Setup initial conditions
#-------------------------------------------------------------------------------
from math import tan
depth = 4.6 #depth of water domain; needs correction is actually slightly <4.6
b = 9.144 #from Synolakis paper
def topography(x,y):#sloping intially (gradient=x/2when b<9.144) then flat (depth=4.6m)
    z = []
    for xi in x:
        if xi < b: 
            z.append(-xi/2.)
        else:
            z.append(-depth)
        #print 'hello', xi,b, z
    return z
domain.set_quantity('stage', 0.0)
domain.set_quantity('friction', 0.0) 
domain.set_quantity('elevation', topography,
                    use_cache = True,
                    verbose = True,
                    alpha = 0.1 #inclination of domain from strike plane
                    )

#-------------------------------------------------------------------------------
# Set up scenario (tsunami_source is a callable object used with set_quantity)
#-------------------------------------------------------------------------------
from smf import slide_tsunami
from math import atan, degrees

a = .4572 # slide height from Synolakis paper
tsunami_source = slide_tsunami(length=0.91,#slide length
                               width=0.61,#slide width
                               depth=.01, # depth below surface to top of slide (delta from Synolakis paper)
                               slope=degrees(atan(0.5)),
                               thickness=a, #slide thickness
                               x0=0.75, 
                               y0=width/2, 
                               alpha=0.,
         		       #dx=0.001,
                               domain=domain,
                               verbose=True)

#-------------------------------------------------------------------------------                                 
# Setup boundary conditions
#-------------------------------------------------------------------------------
print 'Available boundary tags', domain.get_boundary_tags()

Br = Reflective_boundary(domain)
Bd = Dirichlet_boundary([0,0,0])

domain.set_boundary( {'top': Bd,  'bottom': Bd, 'right': Bd, 'left': Bd} )


#-------------------------------------------------------------------------------                                 
# Evolve system through time
#-------------------------------------------------------------------------------
import time
from Numeric import allclose
from anuga.abstract_2d_finite_volumes.quantity import Quantity

t0 = time.time()
finaltime = 500
for t in domain.evolve(yieldstep = 1, finaltime = project.finaltime): 
    domain.write_time()
    stagestep = domain.get_quantity('stage') 

    if allclose(t, 1):
        slide = Quantity(domain)
        slide.set_values(tsunami_source)
        domain.set_quantity('stage', slide + stagestep)
    
print 'That took %.2f seconds' %(time.time()-t0)

print 'finished'