source: trunk/anuga_work/development/Rudy_vandrie_2007/kitchen_sink_1.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 Numeric import array
#from math import sin, pi

#------------------------------------------------------------------------------
# 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.015)    # Constant stage
domain.set_quantity('friction', 0.005) # Constant friction 


#------------------------------------------------------------------------------
# Setup specialised forcing terms
#------------------------------------------------------------------------------
    
def tap(domain):
    """Implement forcing function for "water tap" in the center
    """
    
    flow = min(4*domain.get_time(), 5) # Turn tap up gradually
    stage_update = domain.quantities['stage'].explicit_update

    # Determine indices in tap area
    if not hasattr(domain,'tap_indices'):
        N = len(domain)    
        domain.tap_indices = []
        coordinates = domain.get_centroid_coordinates() 
        for k in range(N):
            x, y = coordinates[k,:] # Centroid
            if ((x-0.5)**2+(y-0.5)**2) < 0.001:     
                domain.tap_indices.append(k)    

    # Update inflow
    for k in domain.tap_indices:
        stage_update[k] += flow                 


def sink(domain):
    """ Implement forcing function for "water sink"
    """
    
    flow = min(2*domain.get_time(), 0.2) # Turn drain up gradually
    stage_update = domain.quantities['stage'].explicit_update    

    # Determine indices in sink area
    if not hasattr(domain,'sink_indices'):
        N = len(domain)    
        domain.sink_indices = []
        coordinates = domain.get_centroid_coordinates() 
        for k in range(N):
            x, y = coordinates[k,:] # Centroid
            if ((x-0.7)**2+(y-0.4)**2) < 0.005:     
                domain.sink_indices.append(k)    
                
    # Update outflow
    for k in domain.sink_indices:
        stage_update[k] -= flow                 

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
#------------------------------------------------------------------------------
for t in domain.evolve(yieldstep = 0.01, finaltime = 14):
    domain.write_time()