source: trunk/anuga_work/development/2010-projects/anuga_1d/sqpipe/sqpipe_domain.py @ 8253

"""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
"""

import numpy

from anuga_1d.base.generic_domain import Generic_domain

#Shallow water domain
class Sqpipe_domain(Generic_domain):

    def __init__(self, coordinates, conserved_quantities, forcing_terms = [], boundary = None, state = None, update_state_flag = True, bulk_modulus = 1.0, tagged_elements = None):

        evolved_quantities = ['area', 'discharge', 'elevation', 'height', 'velocity','width','top','stage']
        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)

        self.bulk_modulus = bulk_modulus

        # Allocate space for tracking state (if not provided)
        if (state is None):
            self.state = numpy.zeros(self.number_of_elements, numpy.int)
        else:
            self.state = state            

        self.update_state_flag = update_state_flag

        # Get computational parameters
        from anuga_1d.config import minimum_allowed_height, g, h0, epsilon
        self.minimum_allowed_height = minimum_allowed_height
        self.g = g
        self.h0 = h0
        self.epsilon = epsilon

        #forcing terms
        for f in forcing_terms:
            self.forcing_terms.append(f)

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

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

        self.quantities_to_be_stored = conserved_quantities

        self.__doc__ = 'sqpipe_domain'

        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 correctness    
        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 check_integrity(self):
        if (len(self.state) != self.number_of_elements):
                raise Exception("state has invalid length")
        Generic_domain.check_integrity(self)

    def compute_fluxes(self):        
        # Set initial timestep to a large value
        import sys
        timestep = float(sys.maxint)

        # Update state before computing flux (if necessary)
        if self.update_state_flag:
            self.update_state()

        # Import flux method and call it
        from anuga_1d.sqpipe.sqpipe_dual_fixed_width_a_d_comp_flux import compute_fluxes_ext
        self.flux_timestep = compute_fluxes_ext(timestep,self)        

    def update_state(self):
        h = self.quantities['height'].centroid_values
        t = self.quantities['top'].centroid_values

        new_state = numpy.zeros_like(self.state)

        # Update boundary state
        new_state[0] = update_cell_state(0, [0, 1], h, t, self.state)
        N = self.number_of_elements-1
        new_state[N] = update_cell_state(N, [N-1, N], h, t, self.state)

        # Update interior states
        for i in range(1, N-1):
            new_state[i] = update_cell_state(i, [i-1, i+1], h, t, self.state)

        self.state = new_state
        #self.state = numpy.where(h>=t, 1, 0)

    def distribute_to_vertices_and_edges(self):
        # Should be implemented by the derived class
        raise Exception("Not implemented")
        
    def evolve(self, yieldstep = None, finaltime = None, duration = None,
               skip_initial_step = False):
        """Specialisation of basic evolve method from parent class
        """

        #Call basic machinery from parent class
        for t in Generic_domain.evolve(self, yieldstep, finaltime,duration,
                                       skip_initial_step):

            #Pass control on to outer loop for more specific actions
            yield(t)

# Auxillary methods

# Method to update the state of a cell
# cell is the index of the cell to update
# neighbours is a list of the indices of it's neighbours
# (for boundary cells, this could include the cell itself) 
def update_cell_state(cell, neighbours, h, t, s):
    # If height is bigger than top, then pressurise
    if h[cell] >= t[cell]:
        return 1
    else:
        # If the cell was pressurised, and all it's neighbours
        # are pressurised, it remains pressurised even though
        # height is less than top
        # Otherwise, it's free surface
        if all([s[i] for i in neighbours]):
            return s[cell]
        else:
            return 0

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