source: inundation/ga/storm_surge/pyvolution-1d/test_quantity.py @ 1507

#!/usr/bin/env python

import unittest
from math import sqrt, pi


from quantity import *
#from config import epsilon
from Numeric import allclose, array, zeros, Float

        
class TestCase(unittest.TestCase):
    def setUp(self):
        from domain import Domain
        
        a = 0.0
        b = 1.0
        c = 2.0
        d = 2.5
        e = 3.1
        f = 4.0

        self.points          = [a, b, c, d, e, f]
        self.centroids       = [(a+b)/2,(b+c)/2,(c+d)/2,(d+e)/2,(e+f)/2]
        self.vertex_values   = [[1.0,2.0],[2.0,3.0],[3.0,4.0],[4.5,5],[5.5,5.6]]
        self.centroid_values = [[1.5,      2.5,      3.5,      4.75,   5.55]]

        self.domain1 = Domain(self.points[0:2])
        self.domain5 = Domain(self.points)
        
    def tearDown(self):
        pass
        #print "  Tearing down"


    def test_creation(self):
        
        quantity = Quantity(self.domain5, self.vertex_values)
        assert allclose(quantity.centroid_values, self.centroid_values)

        try:
            quantity = Quantity()
        except:
            pass
        else:
            raise 'Should have raised empty quantity exception'


        try:
            quantity = Quantity([1,2])
        except AssertionError:
            pass
        except:
            raise 'Should have raised "missing domain object" error'        
        

    def test_creation_zeros(self):
       
        quantity = Quantity(self.domain1)
        assert allclose(quantity.vertex_values, [[0.,0.]])


        quantity = Quantity(self.domain5)
        assert allclose(quantity.vertex_values, [[0.,0.], [0.,0.],
                                                 [0.,0.], [0.,0.],
                                                 [0.,0.]])        

        
    def test_interpolation(self):
        quantity = Quantity(self.domain1, [[1,2]])
        assert allclose(quantity.centroid_values, 1.5) #Centroid


    def test_interpolation2(self):
        quantity = Quantity(self.domain5,
                            [[1,2], [5,5], [0,9], [-6, 3], [3,4]])
        assert allclose(quantity.centroid_values, [1.5, 5., 4.5, -1.5, 3.5 ]) #Centroid

    def test_set_values(self):
        quantity = Quantity(self.domain5)


        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
                            location = 'vertices')
        assert allclose(quantity.vertex_values,
                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid

        #Test default
        quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]])
        assert allclose(quantity.vertex_values,
                        [[1,2], [5,5], [0,0], [-6, 3], [-2,4]])        
        assert allclose(quantity.centroid_values, [1.5, 5., 0., -1.5, 1.0]) #Centroid

        #Test centroids
        quantity.set_values([1,2,3,4,5], location = 'centroids')
        assert allclose(quantity.centroid_values, [1., 2., 3., 4., 5.]) #Centroid

        #Test exceptions
        try:
            quantity.set_values([[1,2], [5,5], [0,0], [-6, 3], [-2,4]],
                                location = 'bas kamel tuba')
        except:
            pass


        try:
            quantity.set_values([[1,2], [0,0]])
        except AssertionError:
            pass
        except:
            raise 'should have raised AssertionError'

    def test_set_values_const(self):
        quantity = Quantity(self.domain5)

        quantity.set_values(1.0, location = 'vertices')
        assert allclose(quantity.vertex_values,
                        [[1,1], [1,1], [1,1], [1,1], [1,1]])
        assert allclose(quantity.centroid_values, [1, 1, 1, 1, 1]) #Centroid


        quantity.set_values(2.0, location = 'centroids')
        assert allclose(quantity.centroid_values, [2, 2, 2, 2, 2])


    def test_set_values_func(self):
        quantity = Quantity(self.domain5)

        def f(x):
            return x**2

        quantity.set_values(f, location = 'vertices')
        assert allclose(quantity.vertex_values,
                        [[0,1], [1,4], [4,6.25], [6.25,9.61], [9.61,16]])        
        assert allclose(quantity.centroid_values,
                        [0.5, 2.5, 5.125, 7.93, 12.805 ])
        
        quantity.set_values(f, location = 'centroids')
        assert allclose(quantity.centroid_values,
                        [0.25, 1.5**2, 2.25**2, 2.8**2, 3.55**2])        


    def test_gradient(self):
        
        quantity = Conserved_quantity(self.domain5)

        #Set up for a gradient of (3,0) at mid triangle 
        quantity.set_values([2.0, 4.0, 8.0, 2.0, 5.0],
                            location = 'centroids')
        
        #Gradients
        quantity.compute_gradients()
        N = quantity.gradients.shape[0]
        
        G = quantity.gradients
        G0 = zeros(N, Float)
        Q = quantity.centroid_values
        X = quantity.domain.centroids

        def grad0(x0,x1,x2,q0,q1,q2):
            return ((q0-q1)/(x0-x1)*(x2-x1) - (q2-q1)/(x2-x1)*(x0-x1))/(x2-x0)

        def grad1(x0,x1,q0,q1):
            return  (q1 - q0)/(x1 - x0)

        #Check Gradients
        G0[0] = grad1(X[0],X[1],Q[0],Q[1])
        for i in range(1,4):
            G0[i] = grad0(X[i-1],X[i],X[i+1],Q[i-1],Q[i],Q[i+1])
        G0[4] = grad1(X[3],X[4],Q[3],Q[4])

        assert allclose(G,G0)

        quantity.extrapolate_first_order()
        
        assert allclose(quantity.vertex_values, [[2., 2.],
                                                 [4., 4.],
                                                 [8., 8.],
                                                 [2., 2.],
                                                 [5., 5.]])

        quantity.extrapolate_second_order()

        V = quantity.domain.vertices
        for k in range(quantity.domain.number_of_elements):
            G0[k] = G0[k]*(V[k,1]-V[k,0])*0.5

        assert allclose(quantity.vertex_values, [[2.-G0[0], 2.+G0[0]],
                                                 [4.-G0[1], 4.+G0[1]],
                                                 [8.-G0[2], 8.+G0[2]],
                                                 [2.-G0[3], 2.+G0[3]],
                                                 [5.-G0[4], 5.+G0[4]]])



    def test_second_order_extrapolation2(self):
        quantity = Conserved_quantity(self.domain5)        

        #Set up for a gradient of (2), f(x) = 2x
        quantity.set_values([1., 3., 4.5, 5.6, 7.1],
                            location = 'centroids')
        
        #Gradients
        quantity.compute_gradients()
        G = quantity.gradients
        
        allclose(G, 2.0)

        quantity.extrapolate_second_order()

        V = quantity.domain.vertices
        Q = quantity.vertex_values
        
        assert allclose(Q[:,0], 2.0*V[:,0])
        assert allclose(Q[:,1], 2.0*V[:,1])        





        


##     def test_first_order_extrapolator(self):
##         quantity = Conserved_quantity(self.mesh4)                

##         #Test centroids
##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid

##         #Extrapolate
##         quantity.extrapolate_first_order()

##         #Check vertices but not edge values
##         assert allclose(quantity.vertex_values,
##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])


    def test_second_order_limit_extrapolator(self):
        quantity = Conserved_quantity(self.domain5)                        

        #Set up for a gradient of (3,0) at mid triangle 
        quantity.set_values([2.0, 4.0, 8.0, 2.0, 0.0],
                            location = 'centroids')
        
        quantity.extrapolate_second_order()
        quantity.limit()

        #Assert that central triangle is limited by neighbours
        assert quantity.vertex_values[0,0] >= 2.0
        assert quantity.vertex_values[0,0] <= 4.0
        assert quantity.vertex_values[0,1] >= 2.0
        assert quantity.vertex_values[0,1] <= 4.0
        
        assert quantity.vertex_values[1,0] >= 2.0
        assert quantity.vertex_values[1,0] <= 8.0
        assert quantity.vertex_values[1,1] >= 2.0
        assert quantity.vertex_values[1,1] <= 8.0

        assert quantity.vertex_values[2,0] >= 2.0
        assert quantity.vertex_values[2,0] <= 8.0
        assert quantity.vertex_values[2,1] >= 2.0
        assert quantity.vertex_values[2,1] <= 8.0

        assert quantity.vertex_values[3,0] >= 0.0
        assert quantity.vertex_values[3,0] <= 8.0
        assert quantity.vertex_values[3,1] >= 0.0
        assert quantity.vertex_values[3,1] <= 8.0

        assert quantity.vertex_values[4,0] >= 0.0
        assert quantity.vertex_values[4,0] <= 2.0
        assert quantity.vertex_values[4,1] >= 0.0
        assert quantity.vertex_values[4,1] <= 2.0

        #Assert that quantities are conserved
        from Numeric import sum
        for k in range(quantity.centroid_values.shape[0]):
            assert allclose (quantity.centroid_values[k],
                             sum(quantity.vertex_values[k,:])/2)
        
        
        

    def test_limiter(self):
        quantity = Conserved_quantity(self.domain5)

        #Create a deliberate overshoot (e.g. from gradient computation)
        quantity.set_values([[3,4], [5,5], [0,0], [-6, 3], [-2,4]],
                            location = 'vertices')        

        #Limit
        quantity.limit()

        #Assert that central triangle is limited by neighbours
        assert quantity.vertex_values[0,1] >= quantity.centroid_values[0]
        assert quantity.vertex_values[1,0] <= quantity.centroid_values[1]
        
        assert quantity.vertex_values[1,1] >= quantity.centroid_values[1]

        
        #Assert that quantities are conserved
        from Numeric import sum
        for k in range(quantity.centroid_values.shape[0]):
            assert allclose (quantity.centroid_values[k],
                             sum(quantity.vertex_values[k,:])/2)
        


##     def test_distribute_first_order(self):
##         quantity = Conserved_quantity(self.mesh4)        

##         #Test centroids
##         quantity.set_values([1.,2.,3.,4.], location = 'centroids')
##         assert allclose(quantity.centroid_values, [1, 2, 3, 4]) #Centroid


##         #Extrapolate
##         quantity.extrapolate_first_order()

##         #Interpolate
##         quantity.interpolate_from_vertices_to_edges()        

##         assert allclose(quantity.vertex_values,
##                         [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]])
##         assert allclose(quantity.edge_values, [[1,1,1], [2,2,2],
##                                                [3,3,3], [4, 4, 4]])



    def test_update_explicit(self):
        quantity = Conserved_quantity(self.domain5)

        #Test centroids
        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid

        #Set explicit_update
        quantity.explicit_update = array( [1.,1.,1.,1.,1.] )

        #Update with given timestep
        quantity.update(0.1)

        x = array([1, 2, 3, 4, 5]) + array( [.1,.1,.1,.1,.1] )
        assert allclose( quantity.centroid_values, x)

    def test_update_semi_implicit(self):
        quantity = Conserved_quantity(self.domain5)

        #Test centroids
        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid

        #Set semi implicit update
        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )

        #Update with given timestep
        quantity.update(0.1)

        x = array([1, 2, 3, 4, 5])/array( [.9,.9,.9,.9,.9] )
        assert allclose( quantity.centroid_values, x)

    def test_both_updates(self):
        quantity = Conserved_quantity(self.domain5)

        #Test centroids
        quantity.set_values([1.,2.,3.,4.,5.], location = 'centroids')
        assert allclose(quantity.centroid_values, [1, 2, 3, 4, 5]) #Centroid

        #Set explicit_update
        quantity.explicit_update = array( [4.,3.,2.,1.,0.0] )
        
        #Set semi implicit update
        quantity.semi_implicit_update = array( [1.,1.,1.,1.,1.] )

        #Update with given timestep
        quantity.update(0.1)

        x = array([1, 2, 3, 4, 5]) + array( [.4,.3,.2,.1,0.0] )
        x /= array( [.9,.9,.9,.9,.9] )
        assert allclose( quantity.centroid_values, x)        

        
        
#-------------------------------------------------------------
if __name__ == "__main__":
    suite = unittest.makeSuite(TestCase,'test')
    runner = unittest.TextTestRunner(verbosity=2)
    runner.run(suite)