"""
Script for debugging potential problem with the boundary
condition: transmissive_momentum_set_stage

When the wave reflects back towards the boundary excessive values
for momentum are generated leading to degenerate timesteps.

Ole Nielsen and Duncan Gray, GA - 2008 
"""


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

# Related major packages
from anuga.shallow_water import Domain
from anuga.shallow_water import Reflective_boundary
from anuga.shallow_water import Dirichlet_boundary
from anuga.shallow_water import Time_boundary
from anuga.shallow_water import File_boundary
from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary


from anuga.abstract_2d_finite_volumes.util import file_function
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular_cross
from anuga.pmesh.mesh_interface import create_mesh_from_regions



# Scenario specific imports
import create_mesh


def elevation_function(x,y):
    from Numeric import zeros, size, Float
    slope = 4 ## Bit of a magic Number
    
    z = zeros(size(x), Float)
    for i in range(len(x)):
        if x[i] < slope:
            z[i] = 0.0
        else:
            z[i] = (x[i]-slope)*0.1
    return z 


#-------------------------------------------------------------------------
# Model parameters
#-------------------------------------------------------------------------
friction=0.01

basename = 'debug_boundary'

boundary_file = 'boundary.tms'
#boundary_file = 'boundary_mellow.tms'



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

xright  = 19.0
ybottom = 0
ytop    = 0.45
xleft = 0.0

#Very simple mesh
points, elements, boundary = rectangular_cross(200, 5,
                                               len1=xright,
                                               len2=ytop,
                                               origin = (0.0, 0.0))



#-------------------------------------------------------------------------
# Setup computational domain
#-------------------------------------------------------------------------
#domain = Domain(mesh_filename, use_cache = False, verbose = True)
domain = Domain(points, elements, boundary, verbose = True)
print domain.statistics()

domain.set_name(basename)
domain.set_datadir('.')


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

domain.set_quantity('stage', 0.06)
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(boundary_file, domain, verbose=True)


Bts = Transmissive_Momentum_Set_Stage_boundary(domain, function)
Br = Reflective_boundary(domain)
domain.set_boundary( {'top': Br, 'right': Br, 'bottom':Br, 'left': Bts} )


#-------------------------------------------------------------------------
# Evolve system through time
#-------------------------------------------------------------------------

for t in domain.evolve(0.1, 15):
    print domain.timestepping_statistics(track_speeds=False)
    print domain.boundary_statistics(tags=['left'])
    
print 'finished'