source: anuga_work/development/anuga_1d/test_shallow_water.py @ 5536

#!/usr/bin/env python

import unittest
from math import sqrt, pi


from shallow_water_domain import *
from Numeric import allclose, array, ones, Float


class Test_Shallow_Water(unittest.TestCase):
    def setUp(self):
        self.points1 = [0.0, 1.0, 2.0, 3.0]
        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]


        pass
        
    def tearDown(self):
        pass
        #print "  Tearing down"


    def test_creation(self):


        D1 = Domain(self.points1)

        assert allclose(D1.centroids, [0.5, 1.5, 2.5])


    def test_compute_fluxes(self):
        

        D1 = Domain(self.points1)
        D1.set_quantity('stage',2.0)
        D1.set_boundary({'exterior' : Reflective_boundary(D1)})

        stage_ud, xmom_ud = compute_fluxes_python(D1)
        
        D1.compute_fluxes()

        #print D1.quantities['xmomentum'].explicit_update
        #print compute_fluxes_python(D1)

        assert allclose( D1.quantities['stage'].explicit_update, stage_ud )
        assert allclose( D1.quantities['xmomentum'].explicit_update, xmom_ud )


#==============================================================================

def compute_fluxes_python(domain):
    """Compute all fluxes and the timestep suitable for all volumes
    in domain.

    Compute total flux for each conserved quantity using "flux_function"

    Fluxes across each edge are scaled by edgelengths and summed up
    Resulting flux is then scaled by area and stored in
    explicit_update for each of the three conserved quantities
    stage, xmomentum and ymomentum

    The maximal allowable speed computed by the flux_function for each volume
    is converted to a timestep that must not be exceeded. The minimum of
    those is computed as the next overall timestep.

    Post conditions:
      domain.explicit_update is reset to computed flux values
      domain.timestep is set to the largest step satisfying all volumes.
    """

    import sys
    from Numeric import zeros, Float

    N = domain.number_of_elements

    tmp0 = zeros((N,),Float)
    tmp1 = zeros((N,),Float)

    #Shortcuts
    Stage = domain.quantities['stage']
    Xmom = domain.quantities['xmomentum']
#    Ymom = domain.quantities['ymomentum']
    Bed = domain.quantities['elevation']

    #Arrays
    #stage = Stage.edge_values
    #xmom =  Xmom.edge_values
 #   ymom =  Ymom.edge_values
    #bed =   Bed.edge_values
    
    stage = Stage.vertex_values
    xmom =  Xmom.vertex_values
    bed =   Bed.vertex_values
    
    #print 'stage edge values', stage
    #print 'xmom edge values', xmom
    #print 'bed values', bed

    stage_bdry = Stage.boundary_values
    xmom_bdry =  Xmom.boundary_values
    #print 'stage_bdry',stage_bdry
    #print 'xmom_bdry', xmom_bdry
#    ymom_bdry =  Ymom.boundary_values

#    flux = zeros(3, Float) #Work array for summing up fluxes
    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(3):
        for i in range(2):
            #Quantities inside volume facing neighbour i
            #ql[0] = stage[k, i]
            #ql[1] = xmom[k, 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_edges[k,i]
                m = domain.neighbour_vertices[k,i]
                #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
                qr[0] = stage[n, m]
                qr[1] = xmom[n, m]
                zr = bed[n, m]


            #Outward pointing normal vector
            normal = domain.normals[k, i]
        
            #Flux computation using provided function
            #edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)
            #print 'ql',ql
            #print 'qr',qr
            

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

            #print 'edgeflux', edgeflux

            # THIS IS THE LINE TO DEAL WITH LEFT AND RIGHT FLUXES
            # flux = edgefluxleft - edgefluxright
            flux -= edgeflux #* domain.edgelengths[k,i]
            #Update optimal_timestep
            try:
                #timestep = min(timestep, 0.5*domain.radii[k]/max_speed)
                timestep = min(timestep, domain.cfl*0.5*domain.areas[k]/max_speed)
            except ZeroDivisionError:
                pass

        #Normalise by area and store for when all conserved
        #quantities get updated
        flux /= domain.areas[k]

        #Stage.explicit_update[k] = flux[0]
        tmp0[k] = flux[0]
        tmp1[k] = flux[1]


    return tmp0, tmp1



#-------------------------------------------------------------
if __name__ == "__main__":
    suite = unittest.makeSuite(Test_Shallow_Water, 'test')    
    #suite = unittest.makeSuite(Test_Quantity, 'test_set_values_from_file_using_polygon')

    #suite = unittest.makeSuite(Test_Quantity, 'test_set_vertex_values_using_general_interface_with_subset')
    #print "restricted test"
    #suite = unittest.makeSuite(Test_Quantity,'verbose_test_set_values_from_UTM_pts')
    runner = unittest.TextTestRunner()
    runner.run(suite)