#!/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)