"""Simple water flow example using ANUGA

Benchmark problem from teh third International Workshop on Long-Wave
Runup Models
http://www.cee.cornell.edu/longwave/index.cfm?page=benchmark&problem=3

"""


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

from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
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 Transmissive_boundary


#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------
from anuga.pmesh.mesh_interface import create_mesh_from_regions
bounding_polygon = [[0,10],[10,10],[10,0],[0,0]]
remainder_res = 0.25
meshname = 'test.msh'
create_mesh_from_regions(bounding_polygon,
                         boundary_tags={'top': [0],
                                        'right': [1],
                                        'bottom': [2],
                                        'left': [3]},
                         maximum_triangle_area=remainder_res,
                         filename=meshname,
                         #interior_regions=interior_regions,
                         use_cache=False,
                         verbose=True)

from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance
from anuga.caching import cache

domain = cache(pmesh_to_domain_instance,
               (meshname, Domain),
               dependencies = [meshname]) # Create domain
domain.set_name('slide')                    # Output to file runup.sww
domain.set_datadir('.')                     # Use current directory for output
domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities
                                    'xmomentum',
                                    'ymomentum'])   


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

from math import tan, radians
beta = radians(5.7)
delta = 1.
mu = 0.01
g = 9.81

def topography(x,y):
    return x*tan(beta)                       # H(x) = x tan(beta)

domain.set_quantity('elevation', topography) # Use function for elevation
domain.set_quantity('friction', 0.1)         # Constant friction 
domain.set_quantity('stage', 0.0)            # Constant negative initial stage

from Numeric import sqrt, exp
t = 0
def ho(x,y):
    return delta*exp(x)

    #return delta*( -(2*sqrt(x*mu*mu/(delta*tan(beta))) - sqrt(g/delta)*mu*t)**2)
#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

from math import sin, pi, exp
Br = Reflective_boundary(domain)      # Solid reflective wall
Bt = Transmissive_boundary(domain)    # Continue all values on boundary 
Bd = Dirichlet_boundary([-0.2,0.,0.]) # Constant boundary values
Bw = Time_boundary(domain=domain,     # Time dependent boundary  
                   f=lambda t: [(.1*sin(t*2*pi)-0.3) * exp(-t), 0.0, 0.0])

# Associate boundary tags with boundary objects
domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})


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

for t in domain.evolve(yieldstep = 0.1, finaltime = 10.0):
    t1 = domain.get_time()
    t = t1
    domain.set_quantity('elevation', ho)
    domain.write_time()