source: anuga_work/development/pipeflow/test_sww_domain.py @ 7840

#!/usr/bin/env python

import unittest
from math import sqrt, pi


from sww_domain import *
import sww_python


import numpy


class Test_Shallow_Water(unittest.TestCase):
    def setUp(self):
        self.points = [0.0, 1.0, 2.0, 3.0]
        self.vertex_values = [[1.0,2.0],[4.0,5.0],[-1.0,2.0]]
        self.points2 = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]
        
    def tearDown(self):
        pass
        #print "  Tearing down"


    def test_creation(self):
        domain = Domain(self.points)
        assert numpy.allclose(domain.centroids, [0.5, 1.5, 2.5])

    def test_compute_fluxes(self):
        """
        Compare shallow_water_domain flux calculation against a previous
        Python implementation (defined in this file)
        """
        domain = Domain(self.points)
        domain.set_quantity('stage',2.0)
        domain.set_boundary({'exterior' : Reflective_boundary(domain)})

        stage_ud, xmom_ud = sww_python.compute_fluxes(domain)

        domain.distribute_to_vertices_and_edges()
        domain.compute_fluxes()

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


    def test_sww_python_flux_function(self):
        normal = 1.0
        ql = numpy.array([1.0, 2.0],numpy.float)
        qr = numpy.array([1.0, 2.0],numpy.float)
        zl = 0.0
        zr = 0.0

        #This assumes h0 = 1.0e-3!!
        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)
        #print maxspeed
        #print edgeflux
        
        assert numpy.allclose([2.0, 8.9], edgeflux, rtol=1.0e-005)
        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)

        normal = -1.0
        ql = numpy.array([1.0, 2.0],numpy.float)
        qr = numpy.array([1.0, 2.0],numpy.float)
        zl = 0.0
        zr = 0.0

        edgeflux, maxspeed = sww_python.flux_function(normal, ql,qr,zl,zr)


        #print maxspeed
        #print edgeflux        
        
        assert numpy.allclose([-2.0, -8.9], edgeflux, rtol=1.0e-005)
        assert numpy.allclose(5.1305, maxspeed, rtol=1.0e-005)



    def test_gravity(self):
        """
        Compare shallow_water_domain gravity calculation
        """

        def slope_one(x):
            return x

        domain = Domain(self.points)
        domain.set_quantity('stage',4.0)
        domain.set_quantity('elevation',slope_one)
        domain.set_boundary({'exterior' : Reflective_boundary(domain)})

        gravity(domain)
        
        #print domain.quantities['stage'].vertex_values 
        #print domain.quantities['elevation'].vertex_values 
        #print domain.quantities['xmomentum'].explicit_update       

        assert numpy.allclose( [-34.3, -24.5, -14.7], domain.quantities['xmomentum'].explicit_update )


    def test_evolve_first_order(self):
        """
        Compare still lake solution for various versions of shallow_water_domain
        """
        
        def slope_square(x):
            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)

        domain = Domain(self.points2)
        domain.set_quantity('stage',10.0)
        domain.set_quantity('elevation',slope_square)
        domain.set_boundary({'exterior' : Reflective_boundary(domain)})

        domain.default_order = 1
        domain.set_timestepping_method('euler')
        
        yieldstep=1.0
        finaltime=1.0

        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
            pass
        

        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )


    def test_evolve_euler_second_order_space(self):
        """
        Compare still lake solution for various versions of shallow_water_domain
        """

        def slope_square(x):
            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)

        domain = Domain(self.points2)
        domain.set_quantity('stage',10.0)
        domain.set_quantity('elevation',slope_square)
        domain.set_boundary({'exterior' : Reflective_boundary(domain)})

        domain.default_order = 2
        domain.set_timestepping_method('rk2')
        yieldstep=1.0
        finaltime=1.0

        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
            pass
        
        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )

    def test_evolve_second_order_space_time(self):
        """
        Compare still lake solution for various versions of shallow_water_domain
        """

        def slope_square(x):
            return numpy.maximum(4.0-(x-5.0)*(x-5.0), 0.0)

        domain = Domain(self.points2)
        domain.set_quantity('stage',10.0)
        domain.set_quantity('elevation',slope_square)
        domain.set_boundary({'exterior' : Reflective_boundary(domain)})

        domain.default_order = 2
        domain.set_timestepping_method('rk3')

        yieldstep=1.0
        finaltime=1.0

        for t in domain.evolve(yieldstep=yieldstep, finaltime=finaltime):
            pass
        
        assert numpy.allclose( 10.0*numpy.ones(10), domain.quantities['stage'].centroid_values )
        assert numpy.allclose( numpy.zeros(10), domain.quantities['xmomentum'].centroid_values )



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

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


    runner = unittest.TextTestRunner()
    runner.run(suite)