#!/usr/bin/env python import unittest from math import sqrt, pi from config import g, epsilon from Numeric import allclose, array, zeros, ones, Float from shallow_water import * #Variable windfield implemented using functions def speed(t,x,y): """Large speeds halfway between center and edges Low speeds at center and edges """ from math import sqrt, exp, cos, pi N = len(x) s = 0*x #New array for k in range(N): r = sqrt(x[k]**2 + y[k]**2) factor = exp( -(r-0.15)**2 ) s[k] = 4000 * factor * (cos(t*2*pi/150) + 2) return s def scalar_func(t,x,y): """Function that returns a scalar. Used to test error message when Numeric array is expected """ return 17.7 def angle(t,x,y): """Rotating field """ from math import sqrt, atan, pi N = len(x) a = 0*x #New array for k in range(N): r = sqrt(x[k]**2 + y[k]**2) angle = atan(y[k]/x[k]) if x[k] < 0: angle+=pi #atan in ]-pi/2; pi/2[ #Take normal direction angle -= pi/2 #Ensure positive radians if angle < 0: angle += 2*pi a[k] = angle/pi*180 return a class TestCase(unittest.TestCase): def setUp(self): pass def tearDown(self): pass def test_rotate(self): normal = array([0.0,-1.0]) q = array([1.0,2.0,3.0]) r = rotate(q, normal, direction = 1) assert r[0] == 1 assert r[1] == -3 assert r[2] == 2 w = rotate(r, normal, direction = -1) assert allclose(w, q) def test_flux_zero_case(self): ql = zeros( 3, Float ) qr = zeros( 3, Float ) normal = zeros( 2, Float ) zl = zr = 0. flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [0,0,0]) assert max_speed == 0. def test_flux_constants(self): w = 2.0 normal = array([1.,0]) ql = array([w, 0, 0]) qr = array([w, 0, 0]) zl = zr = 0. h = w - (zl+zr)/2 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [0., 0.5*g*h**2, 0.]) assert max_speed == sqrt(g*h) def test_flux1(self): #Use data from previous version of pyvolution normal = array([1.,0]) ql = array([-0.2, 2, 3]) qr = array([-0.2, 2, 3]) zl = zr = -0.5 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [2.,13.77433333, 20.]) assert allclose(max_speed, 8.38130948661) def test_flux2(self): #Use data from previous version of pyvolution normal = array([0., -1.]) ql = array([-0.075, 2, 3]) qr = array([-0.075, 2, 3]) zl = zr = -0.375 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [-3.,-20.0, -30.441]) assert allclose(max_speed, 11.7146428199) def test_flux3(self): #Use data from previous version of pyvolution normal = array([-sqrt(2)/2, sqrt(2)/2]) ql = array([-0.075, 2, 3]) qr = array([-0.075, 2, 3]) zl = zr = -0.375 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [sqrt(2)/2, 4.40221112, 7.3829019]) assert allclose(max_speed, 4.0716654239) def test_flux4(self): #Use data from previous version of pyvolution normal = array([-sqrt(2)/2, sqrt(2)/2]) ql = array([-0.34319278, 0.10254161, 0.07273855]) qr = array([-0.30683287, 0.1071986, 0.05930515]) zl = zr = -0.375 flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364]) assert allclose(max_speed, 1.31414103233) def test_sw_domain_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]] domain = Domain(points, vertices) domain.check_integrity() for name in ['level', 'xmomentum', 'ymomentum', 'elevation', 'friction']: assert domain.quantities.has_key(name) assert domain.get_conserved_quantities(0, edge=1) == 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): 'Third', (0, 2): 'First', (2, 0): 'Second', (2, 1): 'Second', (3, 1): 'Second', (3, 2): 'Third'} domain = Domain(points, vertices, boundary) domain.check_integrity() domain.set_quantity('level', [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 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]]) D = Dirichlet_boundary([5,2,1]) T = Transmissive_boundary(domain) R = Reflective_boundary(domain) domain.set_boundary( {'First': D, 'Second': T, 'Third': R}) domain.update_boundary() #Level assert domain.quantities['level'].boundary_values[0] == 2.5 assert domain.quantities['level'].boundary_values[0] ==\ domain.get_conserved_quantities(0, edge=0)[0] #Reflective (2.5) 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] #Reflective (-1.5) #Xmomentum assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective assert domain.quantities['xmomentum'].boundary_values[1] == 2. #Dirichlet assert domain.quantities['xmomentum'].boundary_values[2] ==\ domain.get_conserved_quantities(2, edge=0)[1] #Transmissive assert domain.quantities['xmomentum'].boundary_values[3] ==\ domain.get_conserved_quantities(2, edge=1)[1] #Transmissive assert domain.quantities['xmomentum'].boundary_values[4] ==\ domain.get_conserved_quantities(3, edge=1)[1] #Transmissive assert domain.quantities['xmomentum'].boundary_values[5] == -4.0 #Reflective #Ymomentum assert domain.quantities['ymomentum'].boundary_values[0] == -10.0 #Reflective assert domain.quantities['ymomentum'].boundary_values[1] == 1. #Dirichlet assert domain.quantities['ymomentum'].boundary_values[2] == 30. #Transmissive assert domain.quantities['ymomentum'].boundary_values[3] == 30. #Transmissive assert domain.quantities['ymomentum'].boundary_values[4] == 40. #Transmissive assert domain.quantities['ymomentum'].boundary_values[5] == 40. #Reflective def test_boundary_conditionsII(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', (0, 1): 'Internal'} domain = Domain(points, vertices, boundary) domain.check_integrity() domain.set_quantity('level', [[1,2,3], [5,5,5], [0,0,9], [-6, 3, 3]]) 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]]) D = Dirichlet_boundary([5,2,1]) T = Transmissive_boundary(domain) R = Reflective_boundary(domain) domain.set_boundary( {'First': D, 'Second': T, 'Third': R, 'Internal': None}) domain.update_boundary() domain.check_integrity() def test_compute_fluxes0(self): #Do a full triangle and check that fluxes cancel out for #the constant level case 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]] domain = Domain(points, vertices) domain.set_quantity('level', [[2,2,2], [2,2,2], [2,2,2], [2,2,2]]) domain.check_integrity() zl=zr=0. #Assume flat bed #Flux across right edge of volume 1 normal = domain.get_normal(1,0) ql = domain.get_conserved_quantities(vol_id=1, edge=0) qr = domain.get_conserved_quantities(vol_id=2, edge=2) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) #Check that flux seen from other triangles is inverse tmp = qr; qr=ql; ql=tmp normal = domain.get_normal(2,2) flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux + flux0, 0.) #Flux across upper edge of volume 1 normal = domain.get_normal(1,1) ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) #Check that flux seen from other triangles is inverse tmp = qr; qr=ql; ql=tmp normal = domain.get_normal(3,0) flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux + flux1, 0.) #Flux across lower left hypotenuse of volume 1 normal = domain.get_normal(1,2) ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) #Check that flux seen from other triangles is inverse tmp = qr; qr=ql; ql=tmp normal = domain.get_normal(0,1) flux, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux + flux2, 0.) #Scale by edgelengths, add up anc check that total flux is zero e0 = domain.edgelengths[1, 0] e1 = domain.edgelengths[1, 1] e2 = domain.edgelengths[1, 2] assert allclose(e0*flux0+e1*flux1+e2*flux2, 0.) #Now check that compute_flux yields zeros as well domain.compute_fluxes() for name in ['level', 'xmomentum', 'ymomentum']: #print name, domain.quantities[name].explicit_update assert allclose(domain.quantities[name].explicit_update[1], 0) def test_compute_fluxes1(self): #Use values from previous version 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]] domain = Domain(points, vertices) val0 = 2.+2.0/3 val1 = 4.+4.0/3 val2 = 8.+2.0/3 val3 = 2.+8.0/3 domain.set_quantity('level', [[val0, val0, val0], [val1, val1, val1], [val2, val2, val2], [val3, val3, val3]]) domain.check_integrity() zl=zr=0. #Assume flat bed #Flux across right edge of volume 1 normal = domain.get_normal(1,0) ql = domain.get_conserved_quantities(vol_id=1, edge=0) qr = domain.get_conserved_quantities(vol_id=2, edge=2) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux0, [-15.3598804, 253.71111111, 0.]) assert allclose(max_speed, 9.21592824046) #Flux across upper edge of volume 1 normal = domain.get_normal(1,1) ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux1, [2.4098563, 0., 123.04444444]) assert allclose(max_speed, 7.22956891292) #Flux across lower left hypotenuse of volume 1 normal = domain.get_normal(1,2) ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) assert allclose(flux2, [9.63942522, -61.59685738, -61.59685738]) assert allclose(max_speed, 7.22956891292) #Scale, add up and check that compute_fluxes is correct for vol 1 e0 = domain.edgelengths[1, 0] e1 = domain.edgelengths[1, 1] e2 = domain.edgelengths[1, 2] total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1] assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333]) domain.compute_fluxes() #assert allclose(total_flux, domain.explicit_update[1,:]) for i, name in enumerate(['level', 'xmomentum', 'ymomentum']): assert allclose(total_flux[i], domain.quantities[name].explicit_update[1]) #assert allclose(domain.explicit_update, [ # [0., -69.68888889, -69.68888889], # [-0.68218178, -166.6, -35.93333333], # [-111.77316251, 69.68888889, 0.], # [-35.68522449, 0., 69.68888889]]) assert allclose(domain.quantities['level'].explicit_update, [0., -0.68218178, -111.77316251, -35.68522449]) assert allclose(domain.quantities['xmomentum'].explicit_update, [-69.68888889, -166.6, 69.68888889, 0]) assert allclose(domain.quantities['ymomentum'].explicit_update, [-69.68888889, -35.93333333, 0., 69.68888889]) #assert allclose(domain.quantities[name].explicit_update def test_compute_fluxes2(self): #Random values, incl momentum 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]] domain = Domain(points, vertices) val0 = 2.+2.0/3 val1 = 4.+4.0/3 val2 = 8.+2.0/3 val3 = 2.+8.0/3 zl=zr=0 #Assume flat zero bed domain.set_quantity('elevation', zl*ones( (4,3) )) domain.set_quantity('level', [[val0, val0-1, val0-2], [val1, val1+1, val1], [val2, val2-2, val2], [val3-0.5, val3, val3]]) domain.set_quantity('xmomentum', [[1, 2, 3], [3, 4, 5], [1, -1, 0], [0, -2, 2]]) domain.set_quantity('ymomentum', [[1, -1, 0], [0, -3, 2], [0, 1, 0], [-1, 2, 2]]) domain.check_integrity() #Flux across right edge of volume 1 normal = domain.get_normal(1,0) ql = domain.get_conserved_quantities(vol_id=1, edge=0) qr = domain.get_conserved_quantities(vol_id=2, edge=2) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) #Flux across upper edge of volume 1 normal = domain.get_normal(1,1) ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) #Flux across lower left hypotenuse of volume 1 normal = domain.get_normal(1,2) ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) #Scale, add up and check that compute_fluxes is correct for vol 1 e0 = domain.edgelengths[1, 0] e1 = domain.edgelengths[1, 1] e2 = domain.edgelengths[1, 2] total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1] domain.compute_fluxes() for i, name in enumerate(['level', 'xmomentum', 'ymomentum']): assert allclose(total_flux[i], domain.quantities[name].explicit_update[1]) #assert allclose(total_flux, domain.explicit_update[1,:]) def test_compute_fluxes3(self): #Random values, incl momentum 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]] domain = Domain(points, vertices) val0 = 2.+2.0/3 val1 = 4.+4.0/3 val2 = 8.+2.0/3 val3 = 2.+8.0/3 zl=zr=-3.75 #Assume constant bed (must be less than level) domain.set_quantity('elevation', zl*ones( (4,3) )) domain.set_quantity('level', [[val0, val0-1, val0-2], [val1, val1+1, val1], [val2, val2-2, val2], [val3-0.5, val3, val3]]) domain.set_quantity('xmomentum', [[1, 2, 3], [3, 4, 5], [1, -1, 0], [0, -2, 2]]) domain.set_quantity('ymomentum', [[1, -1, 0], [0, -3, 2], [0, 1, 0], [-1, 2, 2]]) domain.check_integrity() #Flux across right edge of volume 1 normal = domain.get_normal(1,0) ql = domain.get_conserved_quantities(vol_id=1, edge=0) qr = domain.get_conserved_quantities(vol_id=2, edge=2) flux0, max_speed = flux_function(normal, ql, qr, zl, zr) #Flux across upper edge of volume 1 normal = domain.get_normal(1,1) ql = domain.get_conserved_quantities(vol_id=1, edge=1) qr = domain.get_conserved_quantities(vol_id=3, edge=0) flux1, max_speed = flux_function(normal, ql, qr, zl, zr) #Flux across lower left hypotenuse of volume 1 normal = domain.get_normal(1,2) ql = domain.get_conserved_quantities(vol_id=1, edge=2) qr = domain.get_conserved_quantities(vol_id=0, edge=1) flux2, max_speed = flux_function(normal, ql, qr, zl, zr) #Scale, add up and check that compute_fluxes is correct for vol 1 e0 = domain.edgelengths[1, 0] e1 = domain.edgelengths[1, 1] e2 = domain.edgelengths[1, 2] total_flux = -(e0*flux0+e1*flux1+e2*flux2)/domain.areas[1] domain.compute_fluxes() for i, name in enumerate(['level', 'xmomentum', 'ymomentum']): assert allclose(total_flux[i], domain.quantities[name].explicit_update[1]) def test_catching_negative_heights(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]] domain = Domain(points, vertices) val0 = 2.+2.0/3 val1 = 4.+4.0/3 val2 = 8.+2.0/3 val3 = 2.+8.0/3 zl=zr=4 #Too large domain.set_quantity('elevation', zl*ones( (4,3) )) domain.set_quantity('level', [[val0, val0-1, val0-2], [val1, val1+1, val1], [val2, val2-2, val2], [val3-0.5, val3, val3]]) #Should fail try: domain.check_integrity() except: pass ##################################################### def test_gravity(self): #Assuming no friction from config import g 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]] domain = Domain(points, vertices) #Set up for a gradient of (3,0) at mid triangle def slope(x, y): return 3*x h = 0.1 def level(x,y): return slope(x,y)+h domain.set_quantity('elevation', slope) domain.set_quantity('level', level) for name in domain.conserved_quantities: assert allclose(domain.quantities[name].explicit_update, 0) assert allclose(domain.quantities[name].semi_implicit_update, 0) domain.compute_forcing_terms() assert allclose(domain.quantities['level'].explicit_update, 0) assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3) assert allclose(domain.quantities['ymomentum'].explicit_update, 0) def test_manning_friction(self): from config import g 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]] domain = Domain(points, vertices) #Set up for a gradient of (3,0) at mid triangle def slope(x, y): return 3*x h = 0.1 def level(x,y): return slope(x,y)+h eta = 0.07 domain.set_quantity('elevation', slope) domain.set_quantity('level', level) domain.set_quantity('friction', eta) for name in domain.conserved_quantities: assert allclose(domain.quantities[name].explicit_update, 0) assert allclose(domain.quantities[name].semi_implicit_update, 0) domain.compute_forcing_terms() assert allclose(domain.quantities['level'].explicit_update, 0) assert allclose(domain.quantities['xmomentum'].explicit_update, -g*h*3) assert allclose(domain.quantities['ymomentum'].explicit_update, 0) assert allclose(domain.quantities['level'].semi_implicit_update, 0) assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 0) assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0) #Create some momentum for friction to work with domain.set_quantity('xmomentum', 1) S = -g * eta**2 / h**(7.0/3) domain.compute_forcing_terms() assert allclose(domain.quantities['level'].semi_implicit_update, 0) assert allclose(domain.quantities['xmomentum'].semi_implicit_update, S) assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 0) #A more complex example domain.quantities['level'].semi_implicit_update[:] = 0.0 domain.quantities['xmomentum'].semi_implicit_update[:] = 0.0 domain.quantities['ymomentum'].semi_implicit_update[:] = 0.0 domain.set_quantity('xmomentum', 3) domain.set_quantity('ymomentum', 4) S = -g * eta**2 * 5 / h**(7.0/3) domain.compute_forcing_terms() assert allclose(domain.quantities['level'].semi_implicit_update, 0) assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S) assert allclose(domain.quantities['ymomentum'].semi_implicit_update, 4*S) def test_constant_wind_stress(self): from config import rho_a, rho_w, eta_w from math import pi, cos, sin, sqrt 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]] domain = Domain(points, vertices) #Flat surface with 1m of water domain.set_quantity('elevation', 0) domain.set_quantity('level', 1.0) domain.set_quantity('friction', 0) Br = Reflective_boundary(domain) domain.set_boundary({'exterior': Br}) #Setup only one forcing term, constant wind stress s = 100 phi = 135 domain.forcing_terms = [] domain.forcing_terms.append( Wind_stress(s, phi) ) domain.compute_forcing_terms() const = eta_w*rho_a/rho_w #Convert to radians phi = phi*pi/180 #Compute velocity vector (u, v) u = s*cos(phi) v = s*sin(phi) #Compute wind stress S = const * sqrt(u**2 + v**2) assert allclose(domain.quantities['level'].explicit_update, 0) assert allclose(domain.quantities['xmomentum'].explicit_update, S*u) assert allclose(domain.quantities['ymomentum'].explicit_update, S*v) def test_variable_wind_stress(self): from config import rho_a, rho_w, eta_w from math import pi, cos, sin, sqrt 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]] domain = Domain(points, vertices) #Flat surface with 1m of water domain.set_quantity('elevation', 0) domain.set_quantity('level', 1.0) domain.set_quantity('friction', 0) Br = Reflective_boundary(domain) domain.set_boundary({'exterior': Br}) domain.time = 5.54 #Take a random time (not zero) #Setup only one forcing term, constant wind stress s = 100 phi = 135 domain.forcing_terms = [] domain.forcing_terms.append( Wind_stress(s = speed, phi = angle) ) domain.compute_forcing_terms() #Compute reference solution const = eta_w*rho_a/rho_w N = domain.number_of_elements xc = domain.get_centroid_coordinates() t = domain.time x = xc[:,0] y = xc[:,1] s_vec = speed(t,x,y) phi_vec = angle(t,x,y) for k in range(N): #Convert to radians phi = phi_vec[k]*pi/180 s = s_vec[k] #Compute velocity vector (u, v) u = s*cos(phi) v = s*sin(phi) #Compute wind stress S = const * sqrt(u**2 + v**2) assert allclose(domain.quantities['level'].explicit_update[k], 0) assert allclose(domain.quantities['xmomentum'].explicit_update[k], S*u) assert allclose(domain.quantities['ymomentum'].explicit_update[k], S*v) def test_wind_stress_error_condition(self): """Test that windstress reacts properly when forcing functions are wrong - e.g. returns a scalar """ from config import rho_a, rho_w, eta_w from math import pi, cos, sin, sqrt 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]] domain = Domain(points, vertices) #Flat surface with 1m of water domain.set_quantity('elevation', 0) domain.set_quantity('level', 1.0) domain.set_quantity('friction', 0) Br = Reflective_boundary(domain) domain.set_boundary({'exterior': Br}) domain.time = 5.54 #Take a random time (not zero) #Setup only one forcing term, bad func domain.forcing_terms = [] try: domain.forcing_terms.append(Wind_stress(s = scalar_func, phi = angle)) except AssertionError: pass else: msg = 'Should have raised exception' raise msg try: domain.forcing_terms.append(Wind_stress(s = speed, phi = scalar_func)) except AssertionError: pass else: msg = 'Should have raised exception' raise msg try: domain.forcing_terms.append(Wind_stress(s = speed, phi = 'xx')) except: pass else: msg = 'Should have raised exception' raise msg ##################################################### def test_first_order_extrapolator_const_z(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]] domain = Domain(points, vertices) val0 = 2.+2.0/3 val1 = 4.+4.0/3 val2 = 8.+2.0/3 val3 = 2.+8.0/3 zl=zr=-3.75 #Assume constant bed (must be less than level) domain.set_quantity('elevation', zl*ones( (4,3) )) domain.set_quantity('level', [[val0, val0-1, val0-2], [val1, val1+1, val1], [val2, val2-2, val2], [val3-0.5, val3, val3]]) domain.order = 1 domain.distribute_to_vertices_and_edges() #Check that centroid values were distributed to vertices C = domain.quantities['level'].centroid_values for i in range(3): assert allclose( domain.quantities['level'].vertex_values[:,i], C) def test_first_order_limiter_variable_z(self): #Check that first order limiter follows bed_slope from Numeric import alltrue, greater_equal from config import epsilon 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]] domain = Domain(points, vertices) val0 = 2.+2.0/3 val1 = 4.+4.0/3 val2 = 8.+2.0/3 val3 = 2.+8.0/3 domain.set_quantity('elevation', [[0,0,0], [6,0,0], [6,6,6], [6,6,6]]) domain.set_quantity('level', [[val0, val0, val0], [val1, val1, val1], [val2, val2, val2], [val3, val3, val3]]) E = domain.quantities['elevation'].vertex_values L = domain.quantities['level'].vertex_values #Check that some levels are not above elevation (within eps) #- so that the limiter has something to work with assert not alltrue(alltrue(greater_equal(L,E-epsilon))) domain.order = 1 domain.distribute_to_vertices_and_edges() #Check that all levels are above elevation (within eps) assert alltrue(alltrue(greater_equal(L,E-epsilon))) ##################################################### def test_distribute_basic(self): #Using test data generated by pyvolution-2 #Assuming no friction and flat bed (0.0) 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]] domain = Domain(points, vertices) val0 = 2. val1 = 4. val2 = 8. val3 = 2. domain.set_quantity('level', [val0, val1, val2, val3], 'centroids') L = domain.quantities['level'].vertex_values #First order domain.order = 1 domain.distribute_to_vertices_and_edges() assert allclose(L[1], val1) #Second order domain.order = 2 domain.distribute_to_vertices_and_edges() assert allclose(L[1], [2.2, 4.9, 4.9]) def test_distribute_away_from_bed(self): #Using test data generated by pyvolution-2 #Assuming no friction and flat bed (0.0) 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]] domain = Domain(points, vertices) L = domain.quantities['level'].vertex_values def level(x,y): return x**2 domain.set_quantity('level', level, 'centroids') a, b = domain.quantities['level'].compute_gradients() assert allclose(a[1], 3.33333334) assert allclose(b[1], 0.0) domain.order = 1 domain.distribute_to_vertices_and_edges() assert allclose(L[1], 1.77777778) domain.order = 2 domain.distribute_to_vertices_and_edges() assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778]) def test_distribute_away_from_bed1(self): #Using test data generated by pyvolution-2 #Assuming no friction and flat bed (0.0) 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]] domain = Domain(points, vertices) L = domain.quantities['level'].vertex_values def level(x,y): return x**4+y**2 domain.set_quantity('level', level, 'centroids') #print domain.quantities['level'].centroid_values a, b = domain.quantities['level'].compute_gradients() assert allclose(a[1], 25.18518519) assert allclose(b[1], 3.33333333) domain.order = 1 domain.distribute_to_vertices_and_edges() assert allclose(L[1], 4.9382716) domain.order = 2 domain.distribute_to_vertices_and_edges() assert allclose(L[1], [1.07160494, 6.46058131, 7.28262855]) def test_distribute_near_bed(self): #Using test data generated by pyvolution-2 #Assuming no friction and flat bed (0.0) 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]] domain = Domain(points, vertices) #Set up for a gradient of (3,0) at mid triangle def slope(x, y): return 10*x h = 0.1 def level(x,y): return slope(x,y)+h domain.set_quantity('elevation', slope) domain.set_quantity('level', level, 'centroids') #print domain.quantities['elevation'].centroid_values #print domain.quantities['level'].centroid_values E = domain.quantities['elevation'].vertex_values L = domain.quantities['level'].vertex_values #print E domain.order = 1 domain.distribute_to_vertices_and_edges() ##assert allclose(L[1], [0.19999999, 20.05, 20.05]) assert allclose(L[1], [0.1, 20.1, 20.1]) domain.order = 2 domain.distribute_to_vertices_and_edges() assert allclose(L[1], [0.1, 20.1, 20.1]) def test_distribute_near_bed1(self): #Using test data generated by pyvolution-2 #Assuming no friction and flat bed (0.0) 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]] domain = Domain(points, vertices) #Set up for a gradient of (3,0) at mid triangle def slope(x, y): return x**4+y**2 h = 0.1 def level(x,y): return slope(x,y)+h domain.set_quantity('elevation', slope) domain.set_quantity('level', level) #print domain.quantities['elevation'].centroid_values #print domain.quantities['level'].centroid_values E = domain.quantities['elevation'].vertex_values L = domain.quantities['level'].vertex_values #print E domain.order = 1 domain.distribute_to_vertices_and_edges() ##assert allclose(L[1], [4.19999999, 16.07142857, 20.02857143]) assert allclose(L[1], [4.1, 16.1, 20.1]) domain.order = 2 domain.distribute_to_vertices_and_edges() assert allclose(L[1], [4.1, 16.1, 20.1]) def test_second_order_distribute_real_data(self): #Using test data generated by pyvolution-2 #Assuming no friction and flat bed (0.0) a = [0.0, 0.0] b = [0.0, 1.0/5] c = [0.0, 2.0/5] d = [1.0/5, 0.0] e = [1.0/5, 1.0/5] f = [1.0/5, 2.0/5] g = [2.0/5, 2.0/5] points = [a, b, c, d, e, f, g] #bae, efb, cbf, feg vertices = [ [1,0,4], [4,5,1], [2,1,5], [5,4,6]] domain = Domain(points, vertices) def slope(x, y): return -x/3 domain.set_quantity('elevation', slope) domain.set_quantity('level', [0.01298164, 0.00365611, 0.01440365, -0.0381856437096], 'centroids') domain.set_quantity('xmomentum', [0.00670439, 0.01263789, 0.00647805, 0.0178180740668], 'centroids') domain.set_quantity('ymomentum', [-7.23510980e-004, -6.30413883e-005, 6.30413883e-005, 0.000200907255866], 'centroids') E = domain.quantities['elevation'].vertex_values L = domain.quantities['level'].vertex_values X = domain.quantities['xmomentum'].vertex_values Y = domain.quantities['ymomentum'].vertex_values #print E domain.order = 2 domain.distribute_to_vertices_and_edges() #print L[1,:] #print X[1,:] #print Y[1,:] assert allclose(L[1,:], [-0.00825735775384, -0.00801881482869, 0.0272445025825]) assert allclose(X[1,:], [0.0143507718962, 0.0142502147066, 0.00931268339717]) assert allclose(Y[1,:], [-0.000117062180693, 7.94434448109e-005, -0.000151505429018]) def test_second_order_flat_bed_onestep(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary,\ Dirichlet_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=2 # Boundary conditions Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([0.1, 0., 0.]) domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) domain.check_integrity() #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05): pass# domain.write_time() #Data from earlier version of pyvolution assert allclose(domain.min_timestep, 0.0396825396825) assert allclose(domain.max_timestep, 0.0396825396825) assert allclose(domain.quantities['level'].centroid_values[:12], [0.00171396, 0.02656103, 0.00241523, 0.02656103, 0.00241523, 0.02656103, 0.00241523, 0.02656103, 0.00241523, 0.02656103, 0.00241523, 0.0272623]) domain.distribute_to_vertices_and_edges() assert allclose(domain.quantities['level'].vertex_values[:12,0], [0.0001714, 0.02656103, 0.00024152, 0.02656103, 0.00024152, 0.02656103, 0.00024152, 0.02656103, 0.00024152, 0.02656103, 0.00024152, 0.0272623]) assert allclose(domain.quantities['level'].vertex_values[:12,1], [0.00315012, 0.02656103, 0.00024152, 0.02656103, 0.00024152, 0.02656103, 0.00024152, 0.02656103, 0.00024152, 0.02656103, 0.00040506, 0.0272623]) assert allclose(domain.quantities['level'].vertex_values[:12,2], [0.00182037, 0.02656103, 0.00676264, 0.02656103, 0.00676264, 0.02656103, 0.00676264, 0.02656103, 0.00676264, 0.02656103, 0.0065991, 0.0272623]) assert allclose(domain.quantities['xmomentum'].centroid_values[:12], [0.00113961, 0.01302432, 0.00148672, 0.01302432, 0.00148672, 0.01302432, 0.00148672, 0.01302432, 0.00148672 , 0.01302432, 0.00148672, 0.01337143]) assert allclose(domain.quantities['ymomentum'].centroid_values[:12], [-2.91240050e-004 , 1.22721531e-004, -1.22721531e-004, 1.22721531e-004 , -1.22721531e-004, 1.22721531e-004, -1.22721531e-004 , 1.22721531e-004, -1.22721531e-004, 1.22721531e-004, -1.22721531e-004, -4.57969873e-005]) def test_second_order_flat_bed_moresteps(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary,\ Dirichlet_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=2 # Boundary conditions Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([0.1, 0., 0.]) domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) domain.check_integrity() #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): pass #Data from earlier version of pyvolution #assert allclose(domain.min_timestep, 0.0396825396825) #assert allclose(domain.max_timestep, 0.0396825396825) #print domain.quantities['level'].centroid_values def test_flatbed_first_order(self): from mesh_factory import rectangular from shallow_water import Domain,\ Reflective_boundary, Dirichlet_boundary from Numeric import array #Create basic mesh N = 8 points, vertices, boundary = rectangular(N, N) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.visualise = False domain.default_order=1 # Boundary conditions Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([0.2,0.,0.]) domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) domain.check_integrity() #Evolution for t in domain.evolve(yieldstep = 0.02, finaltime = 0.5): pass #domain.write_time() #FIXME: These numbers were from version before 25/10 #assert allclose(domain.min_timestep, 0.0140413643926) #assert allclose(domain.max_timestep, 0.0140947355753) for i in range(3): #assert allclose(domain.quantities['level'].edge_values[:4,i], # [0.10730244,0.12337617,0.11200126,0.12605666]) assert allclose(domain.quantities['xmomentum'].edge_values[:4,i], [0.07610894,0.06901572,0.07284461,0.06819712]) #assert allclose(domain.quantities['ymomentum'].edge_values[:4,i], # [-0.0060238, -0.00157404, -0.00309633, -0.0001637]) def test_flatbed_second_order(self): from mesh_factory import rectangular from shallow_water import Domain,\ Reflective_boundary, Dirichlet_boundary from Numeric import array #Create basic mesh N = 8 points, vertices, boundary = rectangular(N, N) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.visualise = False domain.default_order=2 #domain.minimum_allowed_height = 0.0 #Makes it like the 'oldstyle' balance # Boundary conditions Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([0.2,0.,0.]) domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) domain.check_integrity() #Evolution for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03): pass assert allclose(domain.min_timestep, 0.0210448446782) assert allclose(domain.max_timestep, 0.0210448446782) #print domain.quantities['level'].vertex_values[:4,0] #print domain.quantities['xmomentum'].vertex_values[:4,0] #print domain.quantities['ymomentum'].vertex_values[:4,0] #FIXME: These numbers were from version before 25/10 #assert allclose(domain.quantities['level'].vertex_values[:4,0], # [0.00101913,0.05352143,0.00104852,0.05354394]) assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0], [ 0.00064835, 0.03685719, 0.00085073, 0.03687313]) #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0], # [0.00090581,0.03685719,0.00088303,0.03687313]) assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0], [-0.00139463,0.0006156,-0.00060364,0.00061827]) def test_flatbed_second_order_distribute(self): #Use real data from pyvolution 2 #painfully setup and extracted. from mesh_factory import rectangular from shallow_water import Domain,\ Reflective_boundary, Dirichlet_boundary from Numeric import array #Create basic mesh N = 8 points, vertices, boundary = rectangular(N, N) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.visualise = False domain.default_order=domain.order=2 # Boundary conditions Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([0.2,0.,0.]) domain.set_boundary({'left': Bd, 'right': Br, 'top': Br, 'bottom': Br}) domain.check_integrity() for V in [False, True]: if V: #Set centroids as if system had been evolved L = zeros(2*N*N, Float) L[:32] = [7.21205592e-003, 5.35214298e-002, 1.00910824e-002, 5.35439433e-002, 1.00910824e-002, 5.35439433e-002, 1.00910824e-002, 5.35439433e-002, 1.00910824e-002, 5.35439433e-002, 1.00910824e-002, 5.35439433e-002, 1.00910824e-002, 5.35393928e-002, 1.02344264e-002, 5.59605058e-002, 0.00000000e+000, 3.31027800e-004, 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 4.37962142e-005, 0.00000000e+000, 5.57305948e-005] X = zeros(2*N*N, Float) X[:32] = [6.48351607e-003, 3.68571894e-002, 8.50733285e-003, 3.68731327e-002, 8.50733285e-003, 3.68731327e-002, 8.50733285e-003, 3.68731327e-002, 8.50733285e-003, 3.68731327e-002, 8.50733285e-003, 3.68731327e-002, 8.50733285e-003, 3.68693861e-002, 8.65220973e-003, 3.85055387e-002, 0.00000000e+000, 2.86060840e-004, 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 3.58905503e-005, 0.00000000e+000, 4.57662812e-005] Y = zeros(2*N*N, Float) Y[:32]=[-1.39463104e-003, 6.15600298e-004, -6.03637382e-004, 6.18272251e-004, -6.03637382e-004, 6.18272251e-004, -6.03637382e-004, 6.18272251e-004, -6.03637382e-004, 6.18272251e-004, -6.03637382e-004, 6.18272251e-004, -6.03637382e-004, 6.18599320e-004, -6.74622797e-004, -1.48934756e-004, 0.00000000e+000, -5.35079969e-005, 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, -2.57264987e-005, 0.00000000e+000, -2.57635178e-005] domain.set_quantity('level', L, 'centroids') domain.set_quantity('xmomentum', X, 'centroids') domain.set_quantity('ymomentum', Y, 'centroids') domain.check_integrity() else: #Evolution for t in domain.evolve(yieldstep = 0.01, finaltime = 0.03): pass assert allclose(domain.min_timestep, 0.0210448446782) assert allclose(domain.max_timestep, 0.0210448446782) #Centroids were correct but not vertices. #Hence the check of distribute below. assert allclose(domain.quantities['level'].centroid_values[:4], [0.00721206,0.05352143,0.01009108,0.05354394]) assert allclose(domain.quantities['xmomentum'].centroid_values[:4], [0.00648352,0.03685719,0.00850733,0.03687313]) assert allclose(domain.quantities['ymomentum'].centroid_values[:4], [-0.00139463,0.0006156,-0.00060364,0.00061827]) #print 'C17=', domain.quantities['xmomentum'].centroid_values[17] #print 'C19=', domain.quantities['xmomentum'].centroid_values[19] #assert allclose(domain.quantities['xmomentum'].centroid_values[17],0.00028606084) ##print domain.quantities['xmomentum'].centroid_values[17], V ##print if not V: assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0) else: assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.00028606084) import copy XX = copy.copy(domain.quantities['xmomentum'].centroid_values) assert allclose(domain.quantities['xmomentum'].centroid_values, XX) domain.distribute_to_vertices_and_edges() #assert allclose(domain.quantities['xmomentum'].centroid_values, XX) assert allclose(domain.quantities['xmomentum'].centroid_values[17], 0.0) #FIXME: These numbers were from version before 25/10 #assert allclose(domain.quantities['level'].vertex_values[:4,0], # [0.00101913,0.05352143,0.00104852,0.05354394]) assert allclose(domain.quantities['ymomentum'].vertex_values[:4,0], [-0.00139463,0.0006156,-0.00060364,0.00061827]) assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0], [0.00064835,0.03685719,0.00085073,0.03687313]) #NB NO longer relvant: #This was the culprit. First triangles vertex 0 had an #x-momentum of 0.0064835 instead of 0.00090581 and #third triangle had 0.00850733 instead of 0.00088303 #print domain.quantities['xmomentum'].vertex_values[:4,0] #print domain.quantities['xmomentum'].vertex_values[:4,0] #assert allclose(domain.quantities['xmomentum'].vertex_values[:4,0], # [0.00090581,0.03685719,0.00088303,0.03687313]) def test_bedslope_problem_first_order(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=1 #Bed-slope and friction def x_slope(x, y): return -x/3 domain.set_quantity('elevation', x_slope) # Boundary conditions Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #Initial condition domain.set_quantity('level', Constant_height(x_slope, 0.05)) domain.check_integrity() #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05): pass# domain.write_time() assert allclose(domain.min_timestep, 0.050010003001) assert allclose(domain.max_timestep, 0.050010003001) def test_bedslope_problem_first_order_moresteps(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=1 #Bed-slope and friction def x_slope(x, y): return -x/3 domain.set_quantity('elevation', x_slope) # Boundary conditions Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #Initial condition domain.set_quantity('level', Constant_height(x_slope, 0.05)) domain.check_integrity() #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5): pass# domain.write_time() #Data from earlier version of pyvolution #print domain.quantities['level'].centroid_values assert allclose(domain.quantities['level'].centroid_values, [-0.02998628, -0.01520652, -0.03043492, -0.0149132, -0.03004706, -0.01476251, -0.0298215, -0.01467976, -0.02988158, -0.01474662, -0.03036161, -0.01442995, -0.07624583, -0.06297061, -0.07733792, -0.06342237, -0.07695439, -0.06289595, -0.07635559, -0.0626065, -0.07633628, -0.06280072, -0.07739632, -0.06386738, -0.12161738, -0.11028239, -0.1223796, -0.11095953, -0.12189744, -0.11048616, -0.12074535, -0.10987605, -0.12014311, -0.10976691, -0.12096859, -0.11087692, -0.16868259, -0.15868061, -0.16801135, -0.1588003, -0.16674343, -0.15813323, -0.16457595, -0.15693826, -0.16281096, -0.15585154, -0.16283873, -0.15540068, -0.17450362, -0.19919913, -0.18062882, -0.19764131, -0.17783111, -0.19407213, -0.1736915, -0.19053624, -0.17228678, -0.19105634, -0.17920133, -0.1968828, -0.14244395, -0.14604641, -0.14473537, -0.1506107, -0.14510055, -0.14919522, -0.14175896, -0.14560798, -0.13911658, -0.14439383, -0.13924047, -0.14829043]) def test_bedslope_problem_second_order_one_step(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=2 #Bed-slope and friction at vertices (and interpolated elsewhere) def x_slope(x, y): return -x/3 domain.set_quantity('elevation', x_slope) # Boundary conditions Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #Initial condition domain.set_quantity('level', Constant_height(x_slope, 0.05)) domain.check_integrity() assert allclose(domain.quantities['level'].centroid_values, [0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963]) #print domain.quantities['level'].extrapolate_second_order() #domain.distribute_to_vertices_and_edges() #print domain.quantities['level'].vertex_values[:,0] #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.05): #domain.write_time() pass #print domain.quantities['level'].centroid_values assert allclose(domain.quantities['level'].centroid_values, [0.01290985, 0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.02356019, 0.01619096, 0.0268413, -0.04411074, -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.0248011, -0.04186556, -0.02255593, -0.09966629, -0.08035666, -0.09742112, -0.08035666, -0.09742112, -0.08035666, -0.09742112, -0.08035666, -0.09742112, -0.08035666, -0.09742112, -0.07811149, -0.15522185, -0.13591222, -0.15297667, -0.13591222, -0.15297667, -0.13591222, -0.15297667, -0.13591222, -0.15297667, -0.13591222, -0.15297667, -0.13366704, -0.2107774, -0.19146777, -0.20853223, -0.19146777, -0.20853223, -0.19146777, -0.20853223, -0.19146777, -0.20853223, -0.19146777, -0.20853223, -0.1892226, -0.26120669, -0.24776246, -0.25865535, -0.24776246, -0.25865535, -0.24776246, -0.25865535, -0.24776246, -0.25865535, -0.24776246, -0.25865535, -0.24521113]) def test_bedslope_problem_second_order_two_steps(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=2 #Bed-slope and friction at vertices (and interpolated elsewhere) def x_slope(x, y): return -x/3 domain.set_quantity('elevation', x_slope) # Boundary conditions Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #Initial condition domain.set_quantity('level', Constant_height(x_slope, 0.05)) domain.check_integrity() assert allclose(domain.quantities['level'].centroid_values, [0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963]) #print domain.quantities['level'].extrapolate_second_order() #domain.distribute_to_vertices_and_edges() #print domain.quantities['level'].vertex_values[:,0] #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): pass #Data from earlier version of pyvolution ft=0.1 assert allclose(domain.min_timestep, 0.0376895634803) assert allclose(domain.max_timestep, 0.0415635655309) assert allclose(domain.quantities['level'].centroid_values, [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985, 0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248, 0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789, -0.0229465, -0.04184033, -0.02295693, -0.04184013, -0.02295675, -0.04184486, -0.0228168, -0.04028876, -0.02036486, -0.10029444, -0.08170809, -0.09772846, -0.08021704, -0.09760006, -0.08022143, -0.09759984, -0.08022124, -0.09760261, -0.08008893, -0.09603914, -0.07758209, -0.15584152, -0.13723138, -0.15327266, -0.13572906, -0.15314427, -0.13573349, -0.15314405, -0.13573331, -0.15314679, -0.13560104, -0.15158523, -0.13310701, -0.21208605, -0.19283913, -0.20955631, -0.19134189, -0.20942821, -0.19134598, -0.20942799, -0.1913458, -0.20943005, -0.19120952, -0.20781177, -0.18869401, -0.25384082, -0.2463294, -0.25047649, -0.24464654, -0.25031159, -0.24464253, -0.25031112, -0.24464253, -0.25031463, -0.24454764, -0.24885323, -0.24286438]) def test_bedslope_problem_second_order_two_yieldsteps(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=2 #Bed-slope and friction at vertices (and interpolated elsewhere) def x_slope(x, y): return -x/3 domain.set_quantity('elevation', x_slope) # Boundary conditions Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #Initial condition domain.set_quantity('level', Constant_height(x_slope, 0.05)) domain.check_integrity() assert allclose(domain.quantities['level'].centroid_values, [0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963]) #print domain.quantities['level'].extrapolate_second_order() #domain.distribute_to_vertices_and_edges() #print domain.quantities['level'].vertex_values[:,0] #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): #0.05?? #domain.write_time() pass assert allclose(domain.quantities['level'].centroid_values, [0.00855788, 0.01575204, 0.00994606, 0.01717072, 0.01005985, 0.01716362, 0.01005985, 0.01716299, 0.01007098, 0.01736248, 0.01216452, 0.02026776, -0.04462374, -0.02479045, -0.04199789, -0.0229465, -0.04184033, -0.02295693, -0.04184013, -0.02295675, -0.04184486, -0.0228168, -0.04028876, -0.02036486, -0.10029444, -0.08170809, -0.09772846, -0.08021704, -0.09760006, -0.08022143, -0.09759984, -0.08022124, -0.09760261, -0.08008893, -0.09603914, -0.07758209, -0.15584152, -0.13723138, -0.15327266, -0.13572906, -0.15314427, -0.13573349, -0.15314405, -0.13573331, -0.15314679, -0.13560104, -0.15158523, -0.13310701, -0.21208605, -0.19283913, -0.20955631, -0.19134189, -0.20942821, -0.19134598, -0.20942799, -0.1913458, -0.20943005, -0.19120952, -0.20781177, -0.18869401, -0.25384082, -0.2463294, -0.25047649, -0.24464654, -0.25031159, -0.24464253, -0.25031112, -0.24464253, -0.25031463, -0.24454764, -0.24885323, -0.24286438]) def test_bedslope_problem_second_order_more_steps(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(6, 6) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=2 #Bed-slope and friction at vertices (and interpolated elsewhere) def x_slope(x, y): return -x/3 domain.set_quantity('elevation', x_slope) # Boundary conditions Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #Initial condition domain.set_quantity('level', Constant_height(x_slope, 0.05)) domain.check_integrity() assert allclose(domain.quantities['level'].centroid_values, [0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, 0.01296296, 0.03148148, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.04259259, -0.02407407, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.09814815, -0.07962963, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.1537037 , -0.13518519, -0.1537037, -0.13518519, -0.1537037, -0.13518519, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.20925926, -0.19074074, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963, -0.26481481, -0.2462963]) #print domain.quantities['level'].extrapolate_second_order() #domain.distribute_to_vertices_and_edges() #print domain.quantities['level'].vertex_values[:,0] #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.5): pass assert allclose(domain.quantities['level'].centroid_values, [-0.02907028, -0.01475478, -0.02973417, -0.01447186, -0.02932665, -0.01428285, -0.02901975, -0.0141361, -0.02898816, -0.01418135, -0.02961409, -0.01403487, -0.07597998, -0.06252591, -0.07664854, -0.06312532, -0.07638287, -0.06265139, -0.07571145, -0.06235231, -0.0756817, -0.06245309, -0.07652292, -0.06289946, -0.12367464, -0.11088981, -0.12237277, -0.11115338, -0.1218934, -0.1107174, -0.12081485, -0.11000491, -0.12038451, -0.11010335, -0.12102113, -0.11012105, -0.16909116, -0.15831543, -0.16730214, -0.15786249, -0.1665493, -0.15697919, -0.16496618, -0.15559852, -0.16338679, -0.15509088, -0.16364092, -0.15424423, -0.18771107, -0.19903904, -0.18903759, -0.19858437, -0.18701552, -0.19697797, -0.1833593, -0.19505871, -0.1818806, -0.19418042, -0.18586159, -0.19576946, -0.13986873, -0.14170053, -0.14132188, -0.14560674, -0.14095617, -0.14373292, -0.13785933, -0.14033364, -0.13592955, -0.13936356, -0.13596008, -0.14216296]) assert allclose(domain.quantities['xmomentum'].centroid_values, [ 0.00831121, 0.00317948, 0.00731797, 0.00334939, 0.00764717, 0.00348053, 0.00788729, 0.00356522, 0.00780649, 0.00341919, 0.00693525, 0.00310375, 0.02166196, 0.01421475, 0.02017737, 0.01316839, 0.02037015, 0.01368659, 0.02106, 0.01399161, 0.02074514, 0.01354935, 0.01887407, 0.0123113, 0.03775083, 0.02855197, 0.03689337, 0.02759782, 0.03732848, 0.02812072, 0.03872545, 0.02913348, 0.03880939, 0.02803804, 0.03546499, 0.0260039, 0.0632131, 0.04730634, 0.0576324, 0.04592336, 0.05790921, 0.04690514, 0.05986467, 0.04871165, 0.06170068, 0.04811572, 0.05657041, 0.04416292, 0.08489642, 0.07188097, 0.07835261, 0.06843406, 0.07986412, 0.0698247, 0.08201071, 0.07216756, 0.08378418, 0.07273624, 0.080399, 0.06645841, 0.01631548, 0.04691608, 0.0206632, 0.044409, 0.02115518, 0.04560305, 0.02160608, 0.04663725, 0.02174734, 0.04795559, 0.02281427, 0.05667111]) assert allclose(domain.quantities['ymomentum'].centroid_values, [ 1.45876601e-004, -3.24627393e-004, -1.57572719e-004, -2.92790187e-004, -9.90988382e-005, -3.06677335e-004, -1.62493106e-004, -3.71310004e-004, -1.99445058e-004, -3.28493467e-004, 6.68217349e-005, -8.42042805e-006, 5.05093371e-004, -1.42842214e-004, -6.81454718e-005, -5.02084057e-004, -8.50583861e-005, -4.65443981e-004, -1.96406564e-004, -5.88889562e-004, -2.70160173e-004, -5.35485454e-004, 2.60780997e-004, 3.12145471e-005, 5.16189608e-004, 1.07069062e-004, 9.29989252e-005, -3.71211119e-004, 1.16350246e-004, -3.82407830e-004, -1.62077969e-004, -6.30906636e-004, -4.74025708e-004, -6.94463009e-004, 6.15092843e-005, 2.22106820e-004, -6.29589294e-004, 2.43611937e-004, -5.88125094e-004, -6.94293192e-005, -4.17914641e-004, 6.64609019e-005, -7.68334577e-004, -3.40232101e-004, -1.67424308e-003, -7.39485066e-004, -1.59966988e-003, 5.68262838e-005, -1.48470633e-003, -1.84554882e-003, -2.27200099e-003, -1.67506848e-003, -1.95610258e-003, -1.47638801e-003, -1.73779477e-003, -1.85498791e-003, -2.01357843e-003, -2.17675471e-003, -1.65783870e-003, -1.15818681e-003, -1.18663036e-003, -2.94229849e-003, -3.59309018e-003, -5.13496584e-003, -6.17359400e-003, -5.98761937e-003, -6.00540116e-003, -5.01121966e-003, -4.50964850e-003, -3.06319963e-003, 6.08950810e-004, -4.79537921e-004]) def test_temp_play(self): from mesh_factory import rectangular from shallow_water import Domain, Reflective_boundary, Constant_height from Numeric import array #Create basic mesh points, vertices, boundary = rectangular(5, 5) #Create shallow water domain domain = Domain(points, vertices, boundary) domain.smooth = False domain.default_order=2 #Bed-slope and friction at vertices (and interpolated elsewhere) def x_slope(x, y): return -x/3 domain.set_quantity('elevation', x_slope) # Boundary conditions Br = Reflective_boundary(domain) domain.set_boundary({'left': Br, 'right': Br, 'top': Br, 'bottom': Br}) #Initial condition domain.set_quantity('level', Constant_height(x_slope, 0.05)) domain.check_integrity() #Evolution for t in domain.evolve(yieldstep = 0.05, finaltime = 0.1): pass assert allclose(domain.quantities['level'].centroid_values[:4], [0.00206836, 0.01296714, 0.00363415, 0.01438924]) assert allclose(domain.quantities['xmomentum'].centroid_values[:4], [0.01360154, 0.00671133, 0.01264578, 0.00648503]) assert allclose(domain.quantities['ymomentum'].centroid_values[:4], [-1.19201077e-003, -7.23647546e-004, -6.39083123e-005, 6.29815168e-005]) def test_complex_bed(self): #No friction is tested here from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\ Transmissive_boundary, Time_boundary,\ Weir_simple as Weir, Constant_height from mesh_factory import rectangular from Numeric import array N = 12 points, vertices, boundary = rectangular(N, N/2, len1=1.2,len2=0.6, origin=(-0.07, 0)) domain = Domain(points, vertices, boundary) domain.smooth = False domain.visualise = False domain.default_order=2 inflow_stage = 0.1 Z = Weir(inflow_stage) domain.set_quantity('elevation', Z) Br = Reflective_boundary(domain) Bd = Dirichlet_boundary([inflow_stage, 0.0, 0.0]) domain.set_boundary({'left': Bd, 'right': Br, 'bottom': Br, 'top': Br}) domain.set_quantity('level', Constant_height(Z, 0.)) for t in domain.evolve(yieldstep = 0.02, finaltime = 0.2): pass #print domain.quantities['level'].centroid_values #FIXME: These numbers were from version before 25/10 #assert allclose(domain.quantities['level'].centroid_values, # [3.95822638e-002, 5.61022588e-002, 4.66437868e-002, 5.73081011e-002, # 4.72394613e-002, 5.74684939e-002, 4.74309483e-002, 5.77458084e-002, # 4.80628177e-002, 5.85656225e-002, 4.90498542e-002, 6.02609831e-002, # 1.18470315e-002, 1.75136443e-002, 1.18035266e-002, 2.15565695e-002, # 1.31620268e-002, 2.14351640e-002, 1.32351076e-002, 2.15450687e-002, # 1.36414028e-002, 2.24274619e-002, 1.51689511e-002, 2.21789655e-002, # -7.54337535e-003, -6.86362021e-004, -7.74146760e-003, -1.83756530e-003, # -8.16773628e-003, -4.49916813e-004, -8.08202599e-003, -3.91118720e-004, # -8.10292716e-003, -3.88584984e-004, -7.35226124e-003, 2.73985295e-004, # 1.86166683e-001, 8.74070369e-002, 1.86166712e-001, 8.74035875e-002, # 6.11666935e-002, -3.76173225e-002, -6.38333276e-002, -3.76147365e-002, # 6.11666725e-002, 8.73846774e-002, 1.86166697e-001, 8.74171550e-002, # -4.83333333e-002, 1.18333333e-001, -4.83333333e-002, 1.18333333e-001, # -4.83333333e-002, -6.66666667e-003, -1.73333333e-001, -1.31666667e-001, # -1.73333333e-001, -6.66666667e-003, -4.83333333e-002, 1.18333333e-001, # -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001, # -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001, # -2.48333333e-001, -2.31666667e-001, -2.48333333e-001, -2.31666667e-001, # -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001, # -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001, # -4.65000000e-001, -3.65000000e-001, -4.65000000e-001, -3.65000000e-001, # -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001, # -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001, # -5.98333333e-001, -5.81666667e-001, -5.98333333e-001, -5.81666667e-001, # -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001, # -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001, # -6.48333333e-001, -6.31666667e-001, -6.48333333e-001, -6.31666667e-001, # -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001, # -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001, # -5.31666667e-001, -5.98333333e-001, -5.31666667e-001, -5.98333333e-001, # -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001, # -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001, # -4.98333333e-001, -4.81666667e-001, -4.98333333e-001, -4.81666667e-001, # -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001, # -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001, # -5.48333333e-001, -5.31666667e-001, -5.48333333e-001, -5.31666667e-001]) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(TestCase,'test') #suite = unittest.makeSuite(TestCase,'test_boundary_conditionsII') runner = unittest.TextTestRunner() runner.run(suite)