source: trunk/anuga_work/development/sudi/sw_1d/avalanche/A_velocity/shallow_water_domain_avalanche.py @ 7914

"""Class Domain -
1D interval domains for finite-volume computations of
the shallow water wave equation.

This module contains a specialisation of class Domain from module domain.py
consisting of methods specific to the Shallow Water Wave Equation


U_t + E_x = S

where

U = [w, uh]
E = [uh, u^2h + gh^2/2]
S represents source terms forcing the system
(e.g. gravity, friction, wind stress, ...)

and _t, _x, _y denote the derivative with respect to t, x and y respectiely.

The quantities are

symbol    variable name    explanation
x         x                horizontal distance from origin [m]
z         elevation        elevation of bed on which flow is modelled [m]
h         height           water height above z [m]
w         stage            absolute water level, w = z+h [m]
u                          speed in the x direction [m/s]
uh        xmomentum        momentum in the x direction [m^2/s]

eta                        mannings friction coefficient [to appear]
nu                         wind stress coefficient [to appear]

The conserved quantities are w, uh

For details see e.g.
Christopher Zoppou and Stephen Roberts,
Catastrophic Collapse of Water Supply Reservoirs in Urban Areas,
Journal of Hydraulic Engineering, vol. 127, No. 7 July 1999

John Jakeman, Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
Geoscience Australia, 2006
"""


from domain import *
Generic_Domain = Domain #Rename

#Shallow water domain
class Domain(Generic_Domain):

    def __init__(self, coordinates, boundary = None, tagged_elements = None):
        conserved_quantities = ['stage', 'xmomentum'] #['height', 'xmomentum']
        evolved_quantities = ['stage', 'xmomentum', 'elevation', 'height', 'velocity']
        other_quantities = ['friction']
        Generic_Domain.__init__(self, coordinates, boundary,
                                conserved_quantities, evolved_quantities, other_quantities,
                                tagged_elements)
        
        from config import minimum_allowed_height, g, h0
        self.minimum_allowed_height = minimum_allowed_height
        self.g = g
        self.h0 = h0

        self.forcing_terms.append(gravity_F2)  #gravity, gravity_F1, gravity_F2
        #self.forcing_terms.append(manning_friction)

        #Realtime visualisation
        self.visualiser = None
        self.visualise  = False
        self.visualise_color_stage = False
        self.visualise_stage_range = 1.0
        self.visualise_timer = True
        self.visualise_range_z = None
        
        #Stored output
        self.store = True
        self.format = 'sww'
        self.smooth = True

        #Evolve parametrs
        self.cfl = 1.0
        
        #Reduction operation for get_vertex_values
        from util import mean
        self.reduction = mean
        #self.reduction = min  #Looks better near steep slopes

        self.quantities_to_be_stored = ['stage','xmomentum']

        self.__doc__ = 'shallow_water_domain_avalanche'
        self.check_integrity()


    def set_quantities_to_be_stored(self, q):
        """Specify which quantities will be stored in the sww file.

        q must be either:
          - the name of a quantity
          - a list of quantity names
          - None

        In the two first cases, the named quantities will be stored at each
        yieldstep
        (This is in addition to the quantities elevation and friction)  
        If q is None, storage will be switched off altogether.
        """

        if q is None:
            self.quantities_to_be_stored = []            
            self.store = False
            return

        if isinstance(q, basestring):
            q = [q] # Turn argument into a list

        #Check correcness    
        for quantity_name in q:
            msg = 'Quantity %s is not a valid conserved quantity' %quantity_name
            assert quantity_name in self.conserved_quantities, msg 
        
        self.quantities_to_be_stored = q
        

    def initialise_visualiser(self,scale_z=1.0,rect=None):
        #Realtime visualisation
        if self.visualiser is None:
            from realtime_visualisation_new import Visualiser
            self.visualiser = Visualiser(self,scale_z,rect)
            self.visualiser.setup['elevation']=True
            self.visualiser.updating['stage']=True
        self.visualise = True
        if self.visualise_color_stage == True:
            self.visualiser.coloring['stage'] = True
            self.visualiser.qcolor['stage'] = (0.0, 0.0, 0.8)
        print 'initialise visualiser'
        print self.visualiser.setup
        print self.visualiser.updating

    def check_integrity(self):
        Generic_Domain.check_integrity(self)
        #Check that we are solving the shallow water wave equation

        msg = 'First conserved quantity must be "stage"'
        assert self.conserved_quantities[0] == 'stage', msg
        msg = 'Second conserved quantity must be "xmomentum"'
        assert self.conserved_quantities[1] == 'xmomentum', msg
                
        msg = 'First evolved quantity must be "stage"'
        assert self.evolved_quantities[0] == 'stage', msg
        msg = 'Second evolved quantity must be "xmomentum"'
        assert self.evolved_quantities[1] == 'xmomentum', msg
        msg = 'Third evolved quantity must be "elevation"'
        assert self.evolved_quantities[2] == 'elevation', msg
        msg = 'Fourth evolved quantity must be "height"'
        assert self.evolved_quantities[3] == 'height', msg
        msg = 'Fifth evolved quantity must be "velocity"'
        assert self.evolved_quantities[4] == 'velocity', msg

    def extrapolate_second_order_sw(self):
        #Call correct module function
        #(either from this module or C-extension)
        extrapolate_second_order_sw(self)

    def compute_fluxes(self):
        #Call correct module function(either from this module or C-extension)
        compute_fluxes_C_wellbalanced(self)
        #compute_fluxes_C_nonwellbalanced2(self)

    def compute_timestep(self):
        #Call correct module function
        compute_timestep(self)

    def distribute_to_vertices_and_edges(self):
        #Call correct module function(either from this module or C-extension)
        distribute_to_vertices_and_edges_shv(self)
        #distribute_to_vertices_and_edges_shm(self)
        
    def evolve(self, yieldstep = None, finaltime = None, duration = None, skip_initial_step = False):
        #Call basic machinery from parent class
        for t in Generic_Domain.evolve(self, yieldstep, finaltime, duration, skip_initial_step):
            yield(t)

    def initialise_storage(self):
        """Create and initialise self.writer object for storing data.
        Also, save x and bed elevation
        """

        import data_manager

        #Initialise writer
        self.writer = data_manager.get_dataobject(self, mode = 'w')

        #Store vertices and connectivity
        self.writer.store_connectivity()


    def store_timestep(self, name):
        """Store named quantity and time.

        Precondition:
           self.write has been initialised
        """
        self.writer.store_timestep(name)


#=============== End of Shallow Water Domain ===============================

#Rotation of momentum vector
def rotate(q, normal, direction = 1):
    """Rotate the momentum component q (q[1], q[2])
    from x,y coordinates to coordinates based on normal vector.

    If direction is negative the rotation is inverted.

    Input vector is preserved

    This function is specific to the shallow water wave equation
    """

    from Numeric import Float
    from numpy import zeros    

    assert len(q) == 3,\
           'Vector of conserved quantities must have length 3'\
           'for 2D shallow water equation'

    try:
        l = len(normal)
    except:
        raise 'Normal vector must be an Numeric array'

    assert l == 2, 'Normal vector must have 2 components'


    n1 = normal[0]
    n2 = normal[1]

    r = zeros(len(q), Float) #Rotated quantities
    r[0] = q[0]              #First quantity, height, is not rotated

    if direction == -1:
        n2 = -n2


    r[1] =  n1*q[1] + n2*q[2]
    r[2] = -n2*q[1] + n1*q[2]

    return r


def flux_function(normal, ql, qr, zl, zr):
    """Compute fluxes between volumes for the shallow water wave equation
    cast in terms of w = h+z using the 'central scheme' as described in

    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.

    The implemented formula is given in equation (3.15) on page 714

    Conserved quantities w, uh, are stored as elements 0 and 1
    in the numerical vectors ql an qr.

    Bed elevations zl and zr.
    """

    from config import g, epsilon, h0
    from math import sqrt
    from numpy import array

    #Align momentums with x-axis
    q_left = ql
    q_left[1] = q_left[1]*normal
    q_right = qr
    q_right[1] = q_right[1]*normal
        
    z = 0.5*(zl+zr) #Take average of field values

    w_left  = q_left[0]   #w=h+z
    h_left  = w_left-z
    uh_left = q_left[1]

    if h_left < epsilon:
        u_left = 0.0  #Could have been negative
        h_left = 0.0
    else:
        u_left  = uh_left/(h_left +  h0/h_left)


    uh_left = u_left*h_left


    w_right  = q_right[0]  #w=h+z
    h_right  = w_right-z
    uh_right = q_right[1]

    if h_right < epsilon:
        u_right = 0.0  #Could have been negative
        h_right = 0.0
    else:
        u_right  = uh_right/(h_right + h0/h_right)

    uh_right = u_right*h_right
        

    #We have got h and u at vertex, then the following is the calculation of fluxes!
    soundspeed_left  = sqrt(g*h_left)
    soundspeed_right = sqrt(g*h_right)

    #Maximal wave speed
    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
    
    #Minimal wave speed
    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
    
    #Flux computation
    flux_left  = array([u_left*h_left,
                        u_left*uh_left + 0.5*g*h_left*h_left])
    flux_right = array([u_right*h_right,
                        u_right*uh_right + 0.5*g*h_right*h_right])

    denom = s_max-s_min
    if denom == 0.0:
        edgeflux = array([0.0, 0.0])
        max_speed = 0.0
    else:
        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
        edgeflux += s_max*s_min*(q_right-q_left)/denom
        
        edgeflux[1] = edgeflux[1]*normal

        max_speed = max(abs(s_max), abs(s_min))

    return edgeflux, max_speed

#
def compute_timestep(domain):
    import sys
    from Numeric import Float
    from numpy import zeros    

    N = domain.number_of_elements

    #Shortcuts
    Stage = domain.quantities['stage']
    Xmom = domain.quantities['xmomentum']
    Bed = domain.quantities['elevation']
    
    stage = Stage.vertex_values
    xmom =  Xmom.vertex_values
    bed =   Bed.vertex_values

    stage_bdry = Stage.boundary_values
    xmom_bdry =  Xmom.boundary_values

    flux = zeros(2, Float) #Work array for summing up fluxes
    ql = zeros(2, Float)
    qr = zeros(2, Float)

    #Loop
    timestep = float(sys.maxint)
    enter = True
    for k in range(N):

        flux[:] = 0.  #Reset work array
        for i in range(2):
            #Quantities inside volume facing neighbour i
            ql = [stage[k, i], xmom[k, i]]
            zl = bed[k, i]

            #Quantities at neighbour on nearest face
            n = domain.neighbours[k,i]
            if n < 0:
                m = -n-1 #Convert negative flag to index
                qr[0] = stage_bdry[m]
                qr[1] = xmom_bdry[m]
                zr = zl #Extend bed elevation to boundary
            else:
                m = domain.neighbour_vertices[k,i]
                qr[0] = stage[n, m]
                qr[1] =  xmom[n, m]
                zr = bed[n, m]


            #Outward pointing normal vector
            normal = domain.normals[k, i]

            edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)

            #Update optimal_timestep
            try:
                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
            except ZeroDivisionError:
                pass

    domain.timestep = timestep

# Compute flux definition
# ###################################
def compute_fluxes_C_wellbalanced(domain):
    import sys
    from Numeric import Float
    from numpy import zeros    
    
    N = domain.number_of_elements
    timestep = float(sys.maxint)
    epsilon = domain.epsilon
    g = domain.g
    neighbours = domain.neighbours
    neighbour_vertices = domain.neighbour_vertices
    normals = domain.normals
    areas = domain.areas
    
    Stage    = domain.quantities['stage']
    Xmom     = domain.quantities['xmomentum']
    Bed      = domain.quantities['elevation']
    Height   = domain.quantities['height']
    Velocity = domain.quantities['velocity']


    stage_boundary_values = Stage.boundary_values
    xmom_boundary_values  = Xmom.boundary_values
    bed_boundary_values   = Bed.boundary_values
    height_boundary_values= Height.boundary_values
    vel_boundary_values   = Velocity.boundary_values
    
    stage_explicit_update = Stage.explicit_update
    xmom_explicit_update  = Xmom.explicit_update
    bed_explicit_values   = Bed.explicit_update
    height_explicit_values= Height.explicit_update
    vel_explicit_values   = Velocity.explicit_update
    
    max_speed_array = domain.max_speed_array
   
    domain.distribute_to_vertices_and_edges()
    domain.update_boundary()
    
    stage_V = Stage.vertex_values
    xmom_V  = Xmom.vertex_values
    bed_V   = Bed.vertex_values
    height_V= Height.vertex_values
    vel_V   = Velocity.vertex_values

    number_of_elements = len(stage_V)

    from comp_flux_ext_wellbalanced import compute_fluxes_ext_wellbalanced
    domain.flux_timestep = compute_fluxes_ext_wellbalanced(timestep,
                                  epsilon,
                                  g,
                                                           
                                  neighbours,
                                  neighbour_vertices,
                                  normals,
                                  areas,
                                                           
                                  stage_V,
                                  xmom_V,
                                  bed_V,
                                  height_V,
                                  vel_V,
                                                           
                                  stage_boundary_values,
                                  xmom_boundary_values,
                                  bed_boundary_values,
                                  height_boundary_values,
                                  vel_boundary_values,
                                                           
                                  stage_explicit_update,
                                  xmom_explicit_update,
                                  bed_explicit_values,
                                  height_explicit_values,
                                  vel_explicit_values,
                                                           
                                  number_of_elements,
                                  max_speed_array)


def compute_fluxes_C_nonwellbalanced2(domain):
    import sys
    from Numeric import Float
    from numpy import zeros   
    
    N = domain.number_of_elements
    timestep = float(sys.maxint)
    epsilon = domain.epsilon
    g = domain.g
    neighbours = domain.neighbours
    neighbour_vertices = domain.neighbour_vertices
    normals = domain.normals
    areas = domain.areas
    
    Stage    = domain.quantities['stage']
    Xmom     = domain.quantities['xmomentum']
    Bed      = domain.quantities['elevation']
    Height   = domain.quantities['height']
    Velocity = domain.quantities['velocity']


    stage_boundary_values = Stage.boundary_values
    xmom_boundary_values  = Xmom.boundary_values
    bed_boundary_values   = Bed.boundary_values
    height_boundary_values= Height.boundary_values
    vel_boundary_values   = Velocity.boundary_values
    
    stage_explicit_update = Stage.explicit_update
    xmom_explicit_update  = Xmom.explicit_update
    bed_explicit_values   = Bed.explicit_update
    height_explicit_values= Height.explicit_update
    vel_explicit_values   = Velocity.explicit_update
    
    max_speed_array = domain.max_speed_array
   
    domain.distribute_to_vertices_and_edges()
    domain.update_boundary()
    
    stage_V = Stage.vertex_values
    xmom_V  = Xmom.vertex_values
    bed_V   = Bed.vertex_values
    height_V= Height.vertex_values
    vel_V   = Velocity.vertex_values

    number_of_elements = len(stage_V)

    from comp_flux_ext_nonwellbalanced2 import compute_fluxes_ext_nonwellbalanced2
    domain.flux_timestep = compute_fluxes_ext_nonwellbalanced2(timestep,
                                  epsilon,
                                  g,
                                                           
                                  neighbours,
                                  neighbour_vertices,
                                  normals,
                                  areas,
                                                           
                                  stage_V,
                                  xmom_V,
                                  bed_V,
                                  height_V,
                                  vel_V,
                                                           
                                  stage_boundary_values,
                                  xmom_boundary_values,
                                  bed_boundary_values,
                                  height_boundary_values,
                                  vel_boundary_values,
                                                           
                                  stage_explicit_update,
                                  xmom_explicit_update,
                                  bed_explicit_values,
                                  height_explicit_values,
                                  vel_explicit_values,
                                                           
                                  number_of_elements,
                                  max_speed_array)
 

# ###################################
# Module functions for gradient limiting (distribute_to_vertices_and_edges)

def distribute_to_vertices_and_edges_shv(domain):
    """Distribution from centroids to vertices specific to the
    shallow water wave
    equation.

    It will ensure that h (w-z) is always non-negative even in the
    presence of steep bed-slopes by taking a weighted average between shallow
    and deep cases.

    In addition, all conserved quantities get distributed as per either a
    constant (order==1) or a piecewise linear function (order==2).

    FIXME: more explanation about removal of artificial variability etc

    Precondition:
      All quantities defined at centroids and bed elevation defined at
      vertices.

    Postcondition
      Conserved quantities defined at vertices

    """

    #from config import optimised_gradient_limiter

    #Remove very thin layers of water
    #protect_against_infinitesimal_and_negative_heights(domain)  

    import sys
    from Numeric import Float
    from numpy import zeros
    from config import epsilon, h0

    N = domain.number_of_elements

    #Shortcuts
    Stage = domain.quantities['stage']
    Xmom = domain.quantities['xmomentum']
    Bed = domain.quantities['elevation']
    Height = domain.quantities['height']
    Velocity = domain.quantities['velocity']

    #Arrays   
    w_C   = Stage.centroid_values    
    uh_C  = Xmom.centroid_values    
    z_C   = Bed.centroid_values 
    h_C   = Height.centroid_values
    u_C   = Velocity.centroid_values

    for i in range(N):
        h_C[i] = w_C[i] - z_C[i]                                                
        if h_C[i] <= epsilon:
            uh_C[i] = 0.0
            u_C[i] = 0.0
            w_C[i] = z_C[i]
        else:
            u_C[i]  = uh_C[i]/(h_C[i] + h0/h_C[i])

    for name in ['stage', 'height', 'velocity']:
        Q = domain.quantities[name]
        if domain.order == 1:
            Q.extrapolate_first_order()
        elif domain.order == 2:
            Q.extrapolate_second_order()
        else:
            raise 'Unknown order'

    stage_V = domain.quantities['stage'].vertex_values
    bed_V   = domain.quantities['elevation'].vertex_values
    h_V     = domain.quantities['height'].vertex_values
    u_V     = domain.quantities['velocity'].vertex_values               
    xmom_V  = domain.quantities['xmomentum'].vertex_values 

    bed_V[:] = stage_V - h_V
    xmom_V[:] = u_V * h_V
    
    return



def distribute_to_vertices_and_edges_shm(domain):
    # shm stands for STAGE, HEIGHT, MOMENTUM
    """Distribution from centroids to vertices specific to the
    shallow water wave
    equation.

    It will ensure that h (w-z) is always non-negative even in the
    presence of steep bed-slopes by taking a weighted average between shallow
    and deep cases.

    In addition, all conserved quantities get distributed as per either a
    constant (order==1) or a piecewise linear function (order==2).

    FIXME: more explanation about removal of artificial variability etc

    Precondition:
      All quantities defined at centroids and bed elevation defined at
      vertices.

    Postcondition
      Conserved quantities defined at vertices

    """

    #from config import optimised_gradient_limiter

    #Remove very thin layers of water
    #protect_against_infinitesimal_and_negative_heights(domain)  

    import sys
    from Numeric import Float
    from numpy import zeros
    from config import epsilon, h0

    N = domain.number_of_elements

    #Shortcuts
    Stage = domain.quantities['stage']
    Xmom = domain.quantities['xmomentum']
    Bed = domain.quantities['elevation']
    Height = domain.quantities['height']
    Velocity = domain.quantities['velocity']

    #Arrays   
    w_C   = Stage.centroid_values    
    uh_C  = Xmom.centroid_values    
    z_C   = Bed.centroid_values 
    h_C   = Height.centroid_values
    u_C   = Velocity.centroid_values

    for i in range(N):
        h_C[i] = w_C[i] - z_C[i]                                                
        if h_C[i] <= 0.01:  #epsilon :
            #h_C[i] = 0.0
            #w_C[i] = z_C[i]
            uh_C[i] = 0.0
            u_C[i] = 0.0
            
        if h_C[i] < epsilon:
            h_C[i] = 0.0
            w_C[i] = z_C[i]
        else:
            u_C[i]  = uh_C[i]/h_C[i] #+  h0/h_C[i])

    for name in ['stage', 'height', 'xmomentum']:
        Q = domain.quantities[name]
        if domain.order == 1:
            Q.extrapolate_first_order()
        elif domain.order == 2:
            Q.extrapolate_second_order()
        else:
            raise 'Unknown order'

    stage_V = domain.quantities['stage'].vertex_values
    bed_V   = domain.quantities['elevation'].vertex_values
    h_V     = domain.quantities['height'].vertex_values
    u_V     = domain.quantities['velocity'].vertex_values               
    xmom_V  = domain.quantities['xmomentum'].vertex_values 

    bed_V[:] = stage_V - h_V

    for i in range(N):                                          
        if h_V[i] <= 0.0:
            h_V[i] = 0.0
            stage_V[i] = bed_V[i]
            xmom_V[i] = 0.0
    
    u_V[:]  = xmom_V/(h_V + h0/h_V)    
    
    return




    
#
def protect_against_infinitesimal_and_negative_heights(domain):
    """Protect against infinitesimal heights and associated high velocities
    """

    #Shortcuts
    wc = domain.quantities['stage'].centroid_values
    zc = domain.quantities['elevation'].centroid_values
    xmomc = domain.quantities['xmomentum'].centroid_values
    hc = wc - zc  #Water depths at centroids

    zv = domain.quantities['elevation'].vertex_values
    wv = domain.quantities['stage'].vertex_values
    hv = wv-zv
    xmomv = domain.quantities['xmomentum'].vertex_values
    #remove the above two lines and corresponding code below

    #Update
    for k in range(domain.number_of_elements):

        if hc[k] < domain.minimum_allowed_height:
            #Control stage
            if hc[k] < domain.epsilon:
                wc[k] = zc[k] # Contain 'lost mass' error
                wv[k,0] = zv[k,0]
                wv[k,1] = zv[k,1]
                
            xmomc[k] = 0.0
            
    
def h_limiter(domain):
    """Limit slopes for each volume to eliminate artificial variance
    introduced by e.g. second order extrapolator

    limit on h = w-z

    This limiter depends on two quantities (w,z) so it resides within
    this module rather than within quantity.py
    """

    from Numeric import Float
    from numpy import zeros    

    N = domain.number_of_elements
    beta_h = domain.beta_h

    #Shortcuts
    wc = domain.quantities['stage'].centroid_values
    zc = domain.quantities['elevation'].centroid_values
    hc = wc - zc

    wv = domain.quantities['stage'].vertex_values
    zv = domain.quantities['elevation'].vertex_values
    hv = wv-zv

    hvbar = zeros(hv.shape, Float) #h-limited values

    #Find min and max of this and neighbour's centroid values
    hmax = zeros(hc.shape, Float)
    hmin = zeros(hc.shape, Float)

    for k in range(N):
        hmax[k] = hmin[k] = hc[k]
        for i in range(2):   
            n = domain.neighbours[k,i]
            if n >= 0:
                hn = hc[n] #Neighbour's centroid value

                hmin[k] = min(hmin[k], hn)
                hmax[k] = max(hmax[k], hn)


    #Diffences between centroids and maxima/minima
    dhmax = hmax - hc
    dhmin = hmin - hc

    #Deltas between vertex and centroid values
    dh = zeros(hv.shape, Float)
    for i in range(2):
        dh[:,i] = hv[:,i] - hc

    #Phi limiter
    for k in range(N):

        #Find the gradient limiter (phi) across vertices
        phi = 1.0
        for i in range(2):
            r = 1.0
            if (dh[k,i] > 0): r = dhmax[k]/dh[k,i]
            if (dh[k,i] < 0): r = dhmin[k]/dh[k,i]

            phi = min( min(r*beta_h, 1), phi )

        #Then update using phi limiter
        for i in range(2):
            hvbar[k,i] = hc[k] + phi*dh[k,i]

    return hvbar

def balance_deep_and_shallow(domain):
    """Compute linear combination between stage as computed by
    gradient-limiters limiting using w, and stage computed by
    gradient-limiters limiting using h (h-limiter).
    The former takes precedence when heights are large compared to the
    bed slope while the latter takes precedence when heights are
    relatively small.  Anything in between is computed as a balanced
    linear combination in order to avoid numerical disturbances which
    would otherwise appear as a result of hard switching between
    modes.

    The h-limiter is always applied irrespective of the order.
    """

    #Shortcuts
    wc = domain.quantities['stage'].centroid_values
    zc = domain.quantities['elevation'].centroid_values
    hc = wc - zc

    wv = domain.quantities['stage'].vertex_values
    zv = domain.quantities['elevation'].vertex_values
    hv = wv-zv

    #Limit h
    hvbar = h_limiter(domain)

    for k in range(domain.number_of_elements):
        #Compute maximal variation in bed elevation
        #  This quantitiy is
        #    dz = max_i abs(z_i - z_c)
        #  and it is independent of dimension
        #  In the 1d case zc = (z0+z1)/2
        #  In the 2d case zc = (z0+z1+z2)/3

        dz = max(abs(zv[k,0]-zc[k]),
                 abs(zv[k,1]-zc[k]))#,
#                 abs(zv[k,2]-zc[k]))


        hmin = min( hv[k,:] )

        #Create alpha in [0,1], where alpha==0 means using the h-limited
        #stage and alpha==1 means using the w-limited stage as
        #computed by the gradient limiter (both 1st or 2nd order)

        #If hmin > dz/2 then alpha = 1 and the bed will have no effect
        #If hmin < 0 then alpha = 0 reverting to constant height above bed.

        if dz > 0.0:
            alpha = max( min( 2*hmin/dz, 1.0), 0.0 )
        else:
            #Flat bed
            alpha = 1.0

        alpha  = 0.0
        #Let
        #
        #  wvi be the w-limited stage (wvi = zvi + hvi)
        #  wvi- be the h-limited state (wvi- = zvi + hvi-)
        #
        #
        #where i=0,1,2 denotes the vertex ids
        #
        #Weighted balance between w-limited and h-limited stage is
        #
        #  wvi := (1-alpha)*(zvi+hvi-) + alpha*(zvi+hvi)
        #
        #It follows that the updated wvi is
        #  wvi := zvi + (1-alpha)*hvi- + alpha*hvi
        #
        # Momentum is balanced between constant and limited


        #for i in range(3):
        #    wv[k,i] = zv[k,i] + hvbar[k,i]

        #return

        if alpha < 1:

            #for i in range(3):
            for i in range(2):
                wv[k,i] = zv[k,i] + (1.0-alpha)*hvbar[k,i] + alpha*hv[k,i]

            #Momentums at centroids
            xmomc = domain.quantities['xmomentum'].centroid_values


            #Momentums at vertices
            xmomv = domain.quantities['xmomentum'].vertex_values


            # Update momentum as a linear combination of
            # xmomc and ymomc (shallow) and momentum
            # from extrapolator xmomv and ymomv (deep).
            xmomv[k,:] = (1.0-alpha)*xmomc[k] + alpha*xmomv[k,:]





#########################
#Standard forcing terms:
#
def gravity(domain):
    """Apply gravitational pull in the presence of bed slope
    """

    from util import gradient
    from Numeric import Float
    from numpy import zeros, array, sum

    xmom = domain.quantities['xmomentum'].explicit_update
    stage = domain.quantities['stage'].explicit_update

    Stage = domain.quantities['stage']
    Elevation = domain.quantities['elevation']
    h = Stage.vertex_values - Elevation.vertex_values
    b = Elevation.vertex_values
    w = Stage.vertex_values

    x = domain.get_vertex_coordinates()
    g = domain.g

    for k in range(domain.number_of_elements):
        avg_h = sum( h[k,:] )/2

        #Compute bed slope
        x0, x1 = x[k,:]
        b0, b1 = b[k,:]

        w0, w1 = w[k,:]
        wx = gradient(x0, x1, w0, w1)
        bx = gradient(x0, x1, b0, b1)
        
        #Update momentum (explicit update is reset to source values)
        xmom[k] += -g*bx*avg_h


def gravity_F1(domain):
    """Apply gravitational pull in the presence of bed slope
    """

    from util import gradient
    from Numeric import Float
    from numpy import zeros, array, sum    
    from parameters import F1#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    xmom = domain.quantities['xmomentum'].explicit_update
    stage = domain.quantities['stage'].explicit_update

    Stage = domain.quantities['stage']
    Elevation = domain.quantities['elevation']
    h = Stage.vertex_values - Elevation.vertex_values
    b = Elevation.vertex_values
    w = Stage.vertex_values

    x = domain.get_vertex_coordinates()
    g = domain.g

    for k in range(domain.number_of_elements):
        avg_h = sum( h[k,:] )/2

        #Compute bed slope
        x0, x1 = x[k,:]
        b0, b1 = b[k,:]

        w0, w1 = w[k,:]
        wx = gradient(x0, x1, w0, w1)
        bx = gradient(x0, x1, b0, b1)
        
        #Update momentum (explicit update is reset to source values)
        xmom[k] += -g*bx*avg_h + avg_h*F1#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



def gravity_F2(domain):
    """Apply gravitational pull in the presence of bed slope
    """

    from util import gradient
    from Numeric import Float
    from numpy import zeros, array, sum    
    from parameters import F2#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    xmom = domain.quantities['xmomentum'].explicit_update
    stage = domain.quantities['stage'].explicit_update

    Stage = domain.quantities['stage']
    Elevation = domain.quantities['elevation']
    h = Stage.vertex_values - Elevation.vertex_values
    b = Elevation.vertex_values
    w = Stage.vertex_values

    x = domain.get_vertex_coordinates()
    g = domain.g

    for k in range(domain.number_of_elements):
        avg_h = sum( h[k,:] )/2

        #Compute bed slope
        x0, x1 = x[k,:]
        b0, b1 = b[k,:]

        w0, w1 = w[k,:]
        wx = gradient(x0, x1, w0, w1)
        bx = gradient(x0, x1, b0, b1)
        
        #Update momentum (explicit update is reset to source values)
        xmom[k] += -g*bx*avg_h + avg_h*F2#This is an additional friction!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



 
 
def manning_friction(domain):
    """Apply (Manning) friction to water momentum
    """

    from math import sqrt

    w = domain.quantities['stage'].centroid_values
    z = domain.quantities['elevation'].centroid_values
    h = w-z

    uh = domain.quantities['xmomentum'].centroid_values
    eta = domain.quantities['friction'].centroid_values

    xmom_update = domain.quantities['xmomentum'].semi_implicit_update

    N = domain.number_of_elements
    eps = domain.minimum_allowed_height
    g = domain.g

    for k in range(N):
        if eta[k] >= eps:
            if h[k] >= eps:
                S = -g * eta[k]**2 * uh[k]
                S /= h[k]**(7.0/3)

                #Update momentum
                xmom_update[k] += S*uh[k]

def linear_friction(domain):
    """Apply linear friction to water momentum

    Assumes quantity: 'linear_friction' to be present
    """

    from math import sqrt

    w = domain.quantities['stage'].centroid_values
    z = domain.quantities['elevation'].centroid_values
    h = w-z

    uh = domain.quantities['xmomentum'].centroid_values
    tau = domain.quantities['linear_friction'].centroid_values

    xmom_update = domain.quantities['xmomentum'].semi_implicit_update

    N = domain.number_of_elements
    eps = domain.minimum_allowed_height
    g = domain.g #Not necessary? Why was this added?

    for k in range(N):
        if tau[k] >= eps:
            if h[k] >= eps:
                S = -tau[k]/h[k]

                #Update momentum
                xmom_update[k] += S*uh[k]



def check_forcefield(f):
    """Check that f is either
    1: a callable object f(t,x,y), where x and y are vectors
       and that it returns an array or a list of same length
       as x and y
    2: a scalar
    """

    from Numeric import Float
    from numpy import ones, array    


    if callable(f):
        N = 2
        x = ones(2, Float)
        
        try:
            q = f(1.0, x=x)
        except Exception, e:
            msg = 'Function %s could not be executed:\n%s' %(f, e)
            #FIXME: Reconsider this semantics
            raise msg

        try:
            q = array(q).astype(Float)
        except:
            msg = 'Return value from vector function %s could ' %f
            msg += 'not be converted into a Numeric array of floats.\n'
            msg += 'Specified function should return either list or array.'
            raise msg

        #Is this really what we want?
        msg = 'Return vector from function %s ' %f
        msg += 'must have same lenght as input vectors'
        assert len(q) == N, msg

    else:
        try:
            f = float(f)
        except:
            msg = 'Force field %s must be either a scalar' %f
            msg += ' or a vector function'
            raise msg
    return f

class Wind_stress:
    """Apply wind stress to water momentum in terms of
    wind speed [m/s] and wind direction [degrees]
    """

    def __init__(self, *args, **kwargs):
        """Initialise windfield from wind speed s [m/s]
        and wind direction phi [degrees]

        Inputs v and phi can be either scalars or Python functions, e.g.

        W = Wind_stress(10, 178)

        #FIXME - 'normal' degrees are assumed for now, i.e. the
        vector (1,0) has zero degrees.
        We may need to convert from 'compass' degrees later on and also
        map from True north to grid north.

        Arguments can also be Python functions of t,x,y as in

        def speed(t,x,y):
            ...
            return s

        def angle(t,x,y):
            ...
            return phi

        where x and y are vectors.

        and then pass the functions in

        W = Wind_stress(speed, angle)

        The instantiated object W can be appended to the list of
        forcing_terms as in

        Alternatively, one vector valued function for (speed, angle)
        can be applied, providing both quantities simultaneously.
        As in
        W = Wind_stress(F), where returns (speed, angle) for each t.

        domain.forcing_terms.append(W)
        """

        from config import rho_a, rho_w, eta_w
        from Numeric import Float
        from numpy import array        

        if len(args) == 2:
            s = args[0]
            phi = args[1]
        elif len(args) == 1:
            #Assume vector function returning (s, phi)(t,x,y)
            vector_function = args[0]
            s = lambda t,x: vector_function(t,x=x)[0]
            phi = lambda t,x: vector_function(t,x=x)[1]
        else:
           #Assume info is in 2 keyword arguments

           if len(kwargs) == 2:
               s = kwargs['s']
               phi = kwargs['phi']
           else:
               raise 'Assumes two keyword arguments: s=..., phi=....'

        print 'phi', phi
        self.speed = check_forcefield(s)
        self.phi = check_forcefield(phi)

        self.const = eta_w*rho_a/rho_w


    def __call__(self, domain):
        """Evaluate windfield based on values found in domain
        """

        from math import pi, cos, sin, sqrt
        from Numeric import Float, ArrayType
        from numpy import ones

        xmom_update = domain.quantities['xmomentum'].explicit_update

        N = domain.number_of_elements
        t = domain.time

        if callable(self.speed):
            xc = domain.get_centroid_coordinates()
            s_vec = self.speed(t, xc)
        else:
            #Assume s is a scalar

            try:
                s_vec = self.speed * ones(N, Float)
            except:
                msg = 'Speed must be either callable or a scalar: %s' %self.s
                raise msg


        if callable(self.phi):
            xc = domain.get_centroid_coordinates()
            phi_vec = self.phi(t, xc)
        else:
            #Assume phi is a scalar

            try:
                phi_vec = self.phi * ones(N, Float)
            except:
                msg = 'Angle must be either callable or a scalar: %s' %self.phi
                raise msg

        #assign_windfield_values(xmom_update, ymom_update,
        #                        s_vec, phi_vec, self.const)
        assign_windfield_values(xmom_update, s_vec, phi_vec, self.const)


#def assign_windfield_values(xmom_update, ymom_update,
#                            s_vec, phi_vec, const):
def assign_windfield_values(xmom_update, s_vec, phi_vec, const):
    """Python version of assigning wind field to update vectors.
    A c version also exists (for speed)
    """
    from math import pi, cos, sin, sqrt

    N = len(s_vec)
    for k in range(N):
        s = s_vec[k]
        phi = phi_vec[k]

        #Convert to radians
        phi = phi*pi/180

        #Compute velocity vector (u, v)
        u = s*cos(phi)
        v = s*sin(phi)

        #Compute wind stress
        S = const * u
        xmom_update[k] += S*u