source: trunk/anuga_work/development/Rudy_vandrie_2007/kitchen_sink_2.py @ 7924

"""Example of shallow water wave equation.

Specific methods pertaining to the 2D shallow water equation
are imported from shallow_water
for use with the generic finite volume framework

Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the
numerical vector named conserved_quantities.
"""

#------------------------------------------------------------------------------
# Import necessary modules
#------------------------------------------------------------------------------
from anuga.abstract_2d_finite_volumes.mesh_factory import rectangular
from anuga.shallow_water import Domain
from anuga.shallow_water.shallow_water_domain import Reflective_boundary
from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
from anuga.shallow_water.shallow_water_domain import Inflow

#------------------------------------------------------------------------------
# Setup computational domain
#------------------------------------------------------------------------------

#N = 120
N = 40
points, vertices, boundary = rectangular(N, N)
domain = Domain(points, vertices, boundary)   
domain.set_name('kitchen_sink')                 # Output name
print 'Size', len(domain)


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

domain.set_quantity('elevation', 0.0)  # Zero bed elevation
domain.set_quantity('stage', 0.0)    # Constant stage
domain.set_quantity('friction', 0.005) # Constant friction 


#------------------------------------------------------------------------------
# Setup specialised forcing terms
#------------------------------------------------------------------------------

sink = Inflow(center=[0.7, 0.4], radius=0.0707, flow=0.0)               
tap = Inflow(center=(0.5, 0.5), radius=0.0316,
             flow=lambda t: min(0.01*t, 0.0142)) # Tap turning up       
        

domain.forcing_terms.append(tap)
domain.forcing_terms.append(sink)

#------------------------------------------------------------------------------
# Setup boundary conditions
#------------------------------------------------------------------------------

Br = Reflective_boundary(domain)              # Solid reflective wall
domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br})

#------------------------------------------------------------------------------
# Evolve system through time
#------------------------------------------------------------------------------
from math import sin
for t in domain.evolve(yieldstep = 0.01, finaltime = 20):
    domain.write_time()
    
    if domain.get_time() >= 6 and tap.flow != 0.0:
        print 'Turning tap off'
        tap.flow = 0.0
        
    if domain.get_time() >= 3 and sink.flow == 0.0:
        print 'Turning drain on'
        sink.flow = lambda t: -0.04 + sin(t)*0.005