source: anuga_work/development/flow_1d/sww_flow/sww_vel_domain.py @ 7830

"""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 flow_1d.generic_1d.generic_domain import *
from sww_boundary_conditions import *
from sww_forcing_terms import *

#Shallow water domain
class Domain(Generic_domain):

    def __init__(self, coordinates, boundary = None, tagged_elements = None):

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

        #forcing terms not included in 1d domain ?WHy?
        self.forcing_terms.append(gravity)
        #self.forcing_terms.append(manning_friction)
        #print "\nI have Removed forcing terms line 64 1dsw"

        
        #Stored output
        self.store = True
        self.format = 'sww'
        self.smooth = True

        
        #Reduction operation for get_vertex_values
        from flow_1d.utilities.util import mean
        self.reduction = mean
        #self.reduction = min  #Looks better near steep slopes

        self.set_quantities_to_be_stored(['stage','xmomentum'])

        self.__doc__ = 'sww_vel_domain'

        self.check_integrity()



    def 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

        Generic_domain.check_integrity(self)

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

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

#=============== End of Shallow Water Domain ===============================
#-----------------------------------
# Compute fluxes interface
#-----------------------------------
def compute_fluxes(domain):
    """
    Python version of compute fluxes (local_compute_fluxes)
    is available in test_shallow_water_vel_domain.py
    """
 

    from numpy import zeros 
    import sys

        
    timestep = float(sys.maxint)

    stage = domain.quantities['stage']
    xmom = domain.quantities['xmomentum']
    bed = domain.quantities['elevation']


    from sww_vel_comp_flux_ext import compute_fluxes_ext

    domain.flux_timestep = compute_fluxes_ext(timestep,domain,stage,xmom,bed)


#-----------------------------------
# Compute flux definition with vel
#-----------------------------------
def compute_fluxes_vel(domain):
    from Numeric import zeros, Float
    import sys

        
    timestep = float(sys.maxint)

    stage    = domain.quantities['stage']
    xmom     = domain.quantities['xmomentum']
    bed      = domain.quantities['elevation']
    height   = domain.quantities['height']
    velocity = domain.quantities['velocity']


    from sww_vel_comp_flux_ext import compute_fluxes_vel_ext

    domain.flux_timestep = compute_fluxes_vel_ext(timestep,domain,stage,xmom,bed,height,velocity)



#--------------------------------------------------------------------------
def distribute_to_vertices_and_edges_limit_w_u(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 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
        
    #print id(h_C)
##     for i in range(N):
##      h_C[i] = w_C[i] - z_C[i]                                                
##         if h_C[i] <= 1.0e-12:
##          #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
##          h_C[i] = 0.0
##          w_C[i] = z_C[i]
##          #uh_C[i] = 0.0
            
##  #           u_C[i] = 0.0
##  #       else:
##  #           u_C[i] = uh_C[i]/h_C[i]
                
    h0 = 1.0e-12    
    for i in range(N):
        h_C[i] = w_C[i] - z_C[i]                                                
        if h_C[i] < 1.0e-12:
            u_C[i] = 0.0  #Could have been negative
            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])
            u_C[i]  = uh_C[i]/h_C[i]
        
    for name in [ 'velocity', 'stage' ]:
        Q = domain.quantities[name]
        if domain.order == 1:
            Q.extrapolate_first_order()
        elif domain.order == 2:
            Q.extrapolate_second_order()
        else:
            raise 'Unknown order'

    w_V  = domain.quantities['stage'].vertex_values                     
    z_V  = domain.quantities['elevation'].vertex_values 
    h_V  = domain.quantities['height'].vertex_values
    u_V  = domain.quantities['velocity'].vertex_values          
    uh_V = domain.quantities['xmomentum'].vertex_values 

    #print w_V
    #print z_V
    
    h_V[:,:]    = w_V - z_V
    for i in range(len(h_C)):
        for j in range(2):
            if h_V[i,j] < 0.0 :
                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
                dh = h_V[i,j]
                h_V[i,j] = 0.0
                w_V[i,j] = z_V[i,j]
                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
                
    uh_V[:,:] = u_V * h_V

    
    return

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

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

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

    Postcondition
      Conserved quantities defined at vertices

    """

    import sys
    from Numeric import zeros, Float
    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] <= 1.0e-6:
            #print 'h_C[%d]= %15.5e\n' % (i,h_C[i])
            h_C[i] = 0.0
            w_C[i] = z_C[i]
            uh_C[i] = 0.0
            
    for name in [ 'stage', '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'

    w_V  = domain.quantities['stage'].vertex_values                     
    z_V  = domain.quantities['elevation'].vertex_values 
    h_V  = domain.quantities['height'].vertex_values
    u_V  = domain.quantities['velocity'].vertex_values          
    uh_V = domain.quantities['xmomentum'].vertex_values 

    h_V[:,:]    = w_V - z_V

    for i in range(len(h_C)):
        for j in range(2):
            if h_V[i,j] < 0.0 :
                #print 'h_V[%d,%d] = %f \n' % (i,j,h_V[i,j])                 
                dh = h_V[i,j]
                h_V[i,j] = 0.0
                w_V[i,j] = z_V[i,j]
                h_V[i,(j+1)%2] = h_V[i,(j+1)%2] + dh
                w_V[i,(j+1)%2] = w_V[i,(j+1)%2] + dh
                u_V[i,j] = 0.0
            if h_V[i,j] < h0:
                u_V[i,j]
            else:
                u_V[i,j] = uh_V[i,j]/(h_V[i,j] + h0/h_V[i,j])