#!/usr/bin/env python import unittest from math import sqrt from domain import * from config import epsilon from Numeric import allclose, array class TestCase(unittest.TestCase): def setUp(self): pass def tearDown(self): pass def test_simple(self): a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0,0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe, daf, dae vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]] conserved_quantities = ['level', 'xmomentum', 'ymomentum'] other_quantities = ['elevation', 'friction'] domain = Domain(points, vertices, None, conserved_quantities, other_quantities) domain.check_integrity() for name in conserved_quantities + other_quantities: assert domain.quantities.has_key(name) assert domain.get_conserved_quantities(0, edge=1) == 0. def test_conserved_quantities(self): a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0,0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe, daf, dae vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4], [3,0,5], [3,0,4]] domain = Domain(points, vertices, boundary=None, conserved_quantities =\ ['level', 'xmomentum', 'ymomentum']) domain.set_quantity('level', [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3], [0,0,0], [0,0,0]]) domain.set_quantity('xmomentum', [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3], [0,0,0], [0,0,0]]) domain.check_integrity() #Centroids q = domain.get_conserved_quantities(0) assert allclose(q, [2., 2., 0.]) q = domain.get_conserved_quantities(1) assert allclose(q, [5., 5., 0.]) q = domain.get_conserved_quantities(2) assert allclose(q, [3., 3., 0.]) q = domain.get_conserved_quantities(3) assert allclose(q, [0., 0., 0.]) #Edges q = domain.get_conserved_quantities(0, edge=0) assert allclose(q, [2.5, 2.5, 0.]) q = domain.get_conserved_quantities(0, edge=1) assert allclose(q, [2., 2., 0.]) q = domain.get_conserved_quantities(0, edge=2) assert allclose(q, [1.5, 1.5, 0.]) for i in range(3): q = domain.get_conserved_quantities(1, edge=i) assert allclose(q, [5, 5, 0.]) q = domain.get_conserved_quantities(2, edge=0) assert allclose(q, [4.5, 4.5, 0.]) q = domain.get_conserved_quantities(2, edge=1) assert allclose(q, [4.5, 4.5, 0.]) q = domain.get_conserved_quantities(2, edge=2) assert allclose(q, [0., 0., 0.]) q = domain.get_conserved_quantities(3, edge=0) assert allclose(q, [3., 3., 0.]) q = domain.get_conserved_quantities(3, edge=1) assert allclose(q, [-1.5, -1.5, 0.]) q = domain.get_conserved_quantities(3, edge=2) assert allclose(q, [-1.5, -1.5, 0.]) def test_boundary_conditions(self): a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0,0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] boundary = { (0, 0): 'First', (0, 2): 'First', (2, 0): 'Second', (2, 1): 'Second', (3, 1): 'Second', (3, 2): 'Second'} domain = Domain(points, vertices, boundary, conserved_quantities =\ ['level', 'xmomentum', 'ymomentum']) domain.check_integrity() domain.set_quantity('level', [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) domain.set_boundary( {'First': Dirichlet_boundary([5,2,1]), 'Second': Transmissive_boundary(domain)} ) domain.update_boundary() assert domain.quantities['level'].boundary_values[0] == 5. #Dirichlet assert domain.quantities['level'].boundary_values[1] == 5. #Dirichlet assert domain.quantities['level'].boundary_values[2] ==\ domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5) assert domain.quantities['level'].boundary_values[3] ==\ domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5) assert domain.quantities['level'].boundary_values[4] ==\ domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5) assert domain.quantities['level'].boundary_values[5] ==\ domain.get_conserved_quantities(3, edge=2)[0] #Transmissive (-1.5) def test_distribute_first_order(self): """Domain implements a default first order gradient limiter """ a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0,0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] boundary = { (0, 0): 'Third', (0, 2): 'First', (2, 0): 'Second', (2, 1): 'Second', (3, 1): 'Second', (3, 2): 'Third'} domain = Domain(points, vertices, boundary, conserved_quantities =\ ['level', 'xmomentum', 'ymomentum']) domain.check_integrity() domain.set_quantity('level', [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) assert allclose( domain.quantities['level'].centroid_values, [2,5,3,0] ) domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], [3,3,3], [4, 4, 4]]) domain.set_quantity('ymomentum', [[10,10,10], [20,20,20], [30,30,30], [40, 40, 40]]) domain.distribute_to_vertices_and_edges() #First order extrapolation assert allclose( domain.quantities['level'].vertex_values, [[ 2., 2., 2.], [ 5., 5., 5.], [ 3., 3., 3.], [ 0., 0., 0.]]) def test_update_conserved_quantities(self): a = [0.0, 0.0] b = [0.0, 2.0] c = [2.0,0.0] d = [0.0, 4.0] e = [2.0, 2.0] f = [4.0,0.0] points = [a, b, c, d, e, f] #bac, bce, ecf, dbe vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] boundary = { (0, 0): 'Third', (0, 2): 'First', (2, 0): 'Second', (2, 1): 'Second', (3, 1): 'Second', (3, 2): 'Third'} domain = Domain(points, vertices, boundary, conserved_quantities =\ ['level', 'xmomentum', 'ymomentum']) domain.check_integrity() domain.set_quantity('level', [1,2,3,4], 'centroids') domain.set_quantity('xmomentum', [1,2,3,4], 'centroids') domain.set_quantity('ymomentum', [1,2,3,4], 'centroids') #Assign some values to update vectors #Set explicit_update for name in domain.conserved_quantities: domain.quantities[name].explicit_update = array([4.,3.,2.,1.]) domain.quantities[name].semi_implicit_update = array([1.,1.,1.,1.]) #Update with given timestep (assuming no other forcing terms) domain.timestep = 0.1 domain.update_conserved_quantities() x = array([1, 2, 3, 4]) + array( [.4,.3,.2,.1] ) x /= array( [.9,.9,.9,.9] ) for name in domain.conserved_quantities: assert allclose(domain.quantities[name].centroid_values, x) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(TestCase,'test') runner = unittest.TextTestRunner() runner.run(suite)