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

#!/usr/bin/env python

import unittest
from math import sqrt, pi


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

        
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_boundary_allocation(self):
##         quantity = Conserved_quantity(self.mesh4,
##                             [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]])

##         assert quantity.boundary_values.shape[0] == len(self.mesh4.boundary)


    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.mesh4)

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


##         #gradient bewteen t0 and t1 is undefined as det==0
##         assert a[0] == 0.0
##         assert b[0] == 0.0
##         #The others are OK
##         for i in range(1,4):
##             assert a[i] == 3.0
##             assert b[i] == 0.0


##         quantity.extrapolate_second_order()
        
##         assert allclose(quantity.vertex_values, [[2., 2.,  2.],
##                                                  [0., 6.,  6.],
##                                                  [6., 6., 12.],
##                                                  [0., 0.,  6.]])



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

##         #Set up for a gradient of (3,1), f(x) = 3x+y
##         quantity.set_values([2.0+2.0/3, 4.0+4.0/3, 8.0+2.0/3, 2.0+8.0/3],
##                             location = 'centroids')
        
##         #Gradients
##         a, b = quantity.compute_gradients()

##         #gradient bewteen t0 and t1 is undefined as det==0
##         assert a[0] == 0.0
##         assert b[0] == 0.0
##         #The others are OK
##      163.968025   for i in range(1,4):
##             assert allclose(a[i], 3.0)
##             assert allclose(b[i], 1.0)


##         quantity.extrapolate_second_order()

##         assert allclose(quantity.vertex_values[1,0], 2.0)
##         assert allclose(quantity.vertex_values[1,1], 6.0)        
##         assert allclose(quantity.vertex_values[1,2], 8.0)



## #     def test_limiter(self):

## #         initialise_consecutive_datastructure(points=6+4, elements=4)          
        
## #         a = Point (0.0, 0.0)
## #         b = Point (0.0, 2.0)
## #         c = Point (2.0, 0.0)
## #         d = Point (0.0, 4.0)
## #         e = Point (2.0, 2.0)
## #         f = Point (4.0, 0.0)

## #         #Set up for a gradient of (3,1), f(x) = 3x+y
## #         v1 = Volume(b,a,c,array([0.0,0,0]))        
## #         v2 = Volume(b,c,e,array([1.0,0,0]))
## #         v3 = Volume(e,c,f,array([10.0,0,0]))
## #         v4 = Volume(d,b,e,array([0.0,0,0]))

## #         #Setup neighbour structure
## #         domain = Domain([v1,v2,v3,v4])
## #         domain.precompute()
        
## #         #Lets's check first order first, hey
## #    domain.order = 1
## #    domain.limiter = None
## #         distribute_to_vertices_and_edges(domain)
## #         assert allclose(v2.conserved_quantities_vertex0,
## #                         v2.conserved_quantities_centroid)
## #         assert allclose(v2.conserved_quantities_vertex1,
## #                         v2.conserved_quantities_centroid)
## #         assert allclose(v2.conserved_quantities_vertex2,
## #                         v2.conserved_quantities_centroid)


## #         #Gradient of fitted pwl surface
## #         a, b = compute_gradient(v2.id)
        

## #         assert abs(a[0] - 5.0) < epsilon
## #         assert abs(b[0]) < epsilon        
## #         #assert qminr[0] == 0.0
## #         #assert qmaxr[0] == 10.0
        
## #         #And now for the second order stuff        
## #         # - the full second order extrapolation
## #    domain.order = 2        
## #         distribute_to_vertices_and_edges(domain)


## #         qmin = qmax = v2.conserved_quantities_centroid

## #         qmin = minimum(qmin, v1.conserved_quantities_centroid)
## #         qmax = maximum(qmax, v1.conserved_quantities_centroid)

## #         qmin = minimum(qmin, v3.conserved_quantities_centroid)
## #         qmax = maximum(qmax, v3.conserved_quantities_centroid)

## #         qmin = minimum(qmin, v4.conserved_quantities_centroid)
## #         qmax = maximum(qmax, v4.conserved_quantities_centroid)                    
## #         #assert qminr == qmin
## #         #assert qmaxr == qmax

## #         assert v2.conserved_quantities_vertex0 <= qmax
## #         assert v2.conserved_quantities_vertex0 >= qmin
## #         assert v2.conserved_quantities_vertex1 <= qmax
## #         assert v2.conserved_quantities_vertex1 >= qmin
## #         assert v2.conserved_quantities_vertex2 <= qmax
## #         assert v2.conserved_quantities_vertex2 >= qmin                    
163.968025

## #         #Check that volume has been preserved    

## #         q = v2.conserved_quantities_centroid[0]
## #         w = (v2.conserved_quantities_vertex0[0] +
## #              v2.conserved_quantities_vertex1[0] +
## #              v2.conserved_quantities_vertex2[0])/3

## #         assert allclose(q, w)


        


##     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_extrapolator(self):
##         quantity = Conserved_quantity(self.mesh4)                        

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


##         quantity.extrapolate_second_order()
##         quantity.limit()


##         #Assert that central triangle is limited by neighbours
##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[3,1]
        
##         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
##         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
        
##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
##         assert quantity.vertex_values[1,2] >= quantity.vertex_values[3,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,:])/3)
        
        

        

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

##         #Create a deliberate overshoot (e.g. from gradient computation)
##         quantity.set_values([[3,0,3], [2,2,6], [5,3,8], [8,3,5]])


##         #Limit
##         quantity.limit()

##         #Assert that central triangle is limited by neighbours
##         assert quantity.vertex_values[1,0] >= quantity.vertex_values[0,0]
##         assert quantity.vertex_values[1,0] <= quantity.vertex_values[3,1]
        
##         assert quantity.vertex_values[1,1] <= quantity.vertex_values[2,1]
##         assert quantity.vertex_values[1,1] >= quantity.vertex_values[0,2]
        
##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[2,0]
##         assert quantity.vertex_values[1,2] <= quantity.vertex_values[3,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,:])/3)
        


##     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.mesh4)

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

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

##         #Update with given timestep
##         quantity.update(0.1)

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

##     def test_update_semi_implicit(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

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

##         #Update with given timestep
##         quantity.update(0.1)

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

##     def test_both_updates(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

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

##         #Update with given timestep
##         quantity.update(0.1)

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

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