import unittest, os from math import sqrt, pi from shallow_water_1d import * from Numeric import allclose, array, zeros, ones, Float, take from config import g, epsilon #Variable windfield implemented using functions #def speed(t,x,y): def speed(t,x): """Large speeds halfway between center and edges Low speeds at center and edges """ from math import sqrt, exp, cos, pi x = array(x) #y = array(y) N = len(x) s = 0*x #New array for k in range(N): #r = sqrt(x[k]**2 + y[k]**2) r = x[k] factor = exp( -(r-0.15)**2 ) s[k] = 4000 * factor * (cos(t*2*pi/150) + 2) return s #def scalar_func(t,x,y): def scalar_func(t,x): """Function that returns a scalar. Used to test error message when Numeric array is expected """ return 17.7 #def angle(t,x,y): def angle(t,x): """Rotating field """ from math import sqrt, atan, pi x = array(x) #y = array(y) N = len(x) a = 0*x #New array for k in range(N): #r = sqrt(x[k]**2 + y[k]**2) r =x[k] #angle = atan(y[k]/x[k]) angle = atan(1/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_flux_zero_case(self): #ql = zeros( 3, Float ) #qr = zeros( 3, Float ) ql = zeros( 2, Float ) qr = zeros( 2, Float ) #normal = zeros( 2, Float ) zl = zr = 0. flux, max_speed = flux_function(1.0, ql, qr, zl, zr) #assert allclose(flux, [0,0,0]) assert allclose(flux, [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]) ql = array([w, 0]) qr = array([w, 0]) zl = zr = 0. h = w - (zl+zr)/2 #flux, max_speed = flux_function(normal, ql, qr, zl, zr) flux, max_speed = flux_function(1.0, ql, qr, zl, zr) #assert allclose(flux, [0., 0.5*g*h**2, 0.]) assert allclose(flux, [0., 0.5*g*h**2]) 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]) ql = array([-0.2, 2]) qr = array([-0.2, 2]) zl = zr = -0.5 #flux, max_speed = flux_function(normal, ql, qr, zl, zr) flux, max_speed = flux_function(1.0, ql, qr, zl, zr) #assert allclose(flux, [2.,13.77433333, 20.]) assert allclose(flux, [2.,13.77433333]) assert allclose(max_speed, 8.38130948661) def test_flux2(self): #Use data from previous version of pyvolution """ 2D problem normal = array([0., -1.]) ql = array([-0.075, 2, 3]) qr = array([-0.075, 2, 3]) equivlent to solving 1D ql = array([-0.075, -3]) qr = array([-0.075, -3]) with flux not rotated back i.e flux = [-3.,30.441]) """ ql = array([-0.075, -3]) qr = array([-0.075, -3]) zl = zr = -0.375 #flux, max_speed = flux_function(normal, ql, qr, zl, zr) flux, max_speed = flux_function(1.0, ql, qr, zl, zr) #assert allclose(flux, [-3.,-20.0, -30.441]) assert allclose(flux, [-3.,30.441]) assert allclose(max_speed, 11.7146428199) def test_flux3(self): #def test_flux4(self): #Use data from previous version of pyvolution """ change from 3d case with rotation normal = array([-sqrt(2)/2, sqrt(2)/2]) ql = array([-0.34319278, 0.10254161, 0.07273855]) qr = array([-0.30683287, 0.1071986, 0.05930515]) to 1d case ql = [-0.30683287, -0.03386578] qr = [-0.04072675, 0.03883622] """ ql = array([-0.34319278, -0.02107395]) qr = array([-0.30683287, -0.03386578]) zl = zr = -0.375 #flux, max_speed = flux_function(normal, ql, qr, zl, zr) flux, max_speed = flux_function(1.0, ql, qr, zl, zr) #assert allclose(flux, [-0.04072676, -0.07096636, -0.01604364]) assert allclose(flux, [-0.04072676, 0.03883622]) 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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 f=10.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]] #domain = Domain(points, vertices) domain = Domain(points) domain.check_integrity() #for name in ['stage', 'xmomentum', 'ymomentum', for name in ['stage', 'xmomentum', 'elevation', 'friction']: assert domain.quantities.has_key(name) #assert domain.get_conserved_quantities(0, edge=1) == 0. assert domain.get_conserved_quantities(0, vertex=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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 f=10.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'} #there are points -1 segements with edges (0,1) = (left, right) boundary = {(0,0): 'Start', (4,1): 'Finish'} #boundary = {(0,0): 'exterior', (4,1): 'exterior'} #domain = Domain(points, vertices, boundary) domain = Domain(points, boundary) domain.check_integrity() #Sets stages in each interval(element) #domain.set_quantity('stage', [[1,2,3], [5,5,5], # [0,0,9], [-6, 3, 3]]) domain.set_quantity('stage', [[1,2], [5,5], [0,0], [-6,3], [3,0]]) #domain.set_quantity('xmomentum', [[1,1,1], [2,2,2], # [3,3,3], [4, 4, 4]]) domain.set_quantity('xmomentum', [[1,1], [2,2], [3,3], [4, 4], [5,5]]) #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( {'Start': R, 'Finish': R}) #domain.set_boundary( {'exterior': R, 'exterior': R}) domain.update_boundary() #Stage #assert domain.quantities['stage'].boundary_values[0] == 2.5 assert domain.quantities['stage'].boundary_values[0] ==\ domain.get_conserved_quantities(0, vertex=0)[0] #Reflective (2.5) #assert domain.quantities['stage'].boundary_values[1] == 5. #Dirichlet #assert domain.quantities['stage'].boundary_values[2] ==\ # domain.get_conserved_quantities(2, edge=0)[0] #Transmissive (4.5) #assert domain.quantities['stage'].boundary_values[3] ==\ # domain.get_conserved_quantities(2, edge=1)[0] #Transmissive (4.5) #assert domain.quantities['stage'].boundary_values[4] ==\ # domain.get_conserved_quantities(3, edge=1)[0] #Transmissive (-1.5) #assert domain.quantities['stage'].boundary_values[5] ==\ # domain.get_conserved_quantities(3, edge=2)[0] #Reflective (-1.5) #Xmomentum #print 'xmom' print domain.quantities['xmomentum'].boundary_values[0] assert domain.quantities['xmomentum'].boundary_values[0] == 1.0 #Reflective #SHOULDN'T THIS BE -1.0 #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[1] == -5.0 #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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 f=10.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'} boundary = {(0,0): 'Start', (1,0) : 'Internal', (4,1): 'Finish'} #domain = Domain(points, vertices, boundary) domain = Domain(points, boundary) domain.check_integrity() #domain.set_quantity('stage', [[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]]) domain.set_quantity('stage', [[1,2], [5,5], [0,0], [-6, 3], [3, 0]]) domain.set_quantity('xmomentum', [[1,1], [2,2], [3,3], [4, 4], [5, 5]]) #D = Dirichlet_boundary([5,2,1]) #T = Transmissive_boundary(domain) R = Reflective_boundary(domain) domain.set_boundary( {'Start': R, 'Finish': 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 stage case print 'check min time step in compute fluxes is ok, John' #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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 f=10.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) #domain.set_quantity('stage', [[2,2,2], [2,2,2], # [2,2,2], [2,2,2]]) domain.set_quantity('stage', [[2,2], [2,2], [2,2], [2,2], [2,2]]) domain.check_integrity() #assert allclose(domain.neighbours, [[-1,1,-1], [2,3,0], [-1,-1,1],[1,-1,-1]]) assert allclose(domain.neighbours, [[-1,1], [0,2], [1,3],[2,4], [3,-1]]) #assert allclose(domain.neighbour_edges, [[-1,2,-1], [2,0,1], [-1,-1,0],[1,-1,-1]]) #assert allclose(domain.neighbour_edges, [[-1,0], [1,0], [1,0], [1,0], [1,-1]]) assert allclose(domain.neighbour_vertices, [[-1,0], [1,0], [1,0], [1,0], [1,-1]]) 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) #Flux across right edge of element 1 ql = domain.get_conserved_quantities(vol_id=1, vertex=1) qr = domain.get_conserved_quantities(vol_id=2, vertex=0) #print 'qr and ql 1' #print qr #print ql flux0, max_speed = flux_function(1.0, ql, qr, zl, zr) #Check that flux seen from other triangles is inverse tmp = qr; qr=ql; ql=tmp #print 'qr and ql 2' #print qr #print ql #normal = domain.get_normal(2,2) #flux, max_speed = flux_function(normal, ql, qr, zl, zr) flux, max_speed = flux_function(-1.0, ql, qr, zl, zr) #print 'fluxes' #print flux0 #print flux 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) #Flux across left edge of element 1 #ql = domain.get_conserved_quantities(vol_id=1, edge=0) #qr = domain.get_conserved_quantities(vol_id=0, edge=1) ql = domain.get_conserved_quantities(vol_id=1, vertex=0) qr = domain.get_conserved_quantities(vol_id=0, vertex=1) flux1, max_speed = flux_function(-1.0, 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) flux, max_speed = flux_function(1.0, 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.) print 'flux0',flux0 print 'flux1',flux1 assert allclose(flux0+flux1, 0.) #Now check that compute_flux yields zeros as well domain.compute_fluxes() #for name in ['stage', 'xmomentum', 'ymomentum']: for name in ['stage', 'xmomentum']: #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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) val0 = 2.+2.0/3 val1 = 4.+4.0/3 val2 = 8.+2.0/3 val3 = 2.+8.0/3 #domain.set_quantity('stage', [[val0, val0, val0], [val1, val1, val1], # [val2, val2, val2], [val3, val3, val3]]) domain.set_quantity('stage', [[val0, val0], [val1, val1], [val2, val2], [val3, val3]]) stage = domain.get_quantity('stage') print stage 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) #Flux across right edge of volume 1 #ql = domain.get_conserved_quantities(vol_id=1, edge=1) #qr = domain.get_conserved_quantities(vol_id=2, edge=0) ql = domain.get_conserved_quantities(vol_id=1, vertex=1) qr = domain.get_conserved_quantities(vol_id=2, vertex=0) print 'ql',ql print 'qr',qr flux0, max_speed = flux_function(1.0, ql, qr, zl, zr) print 'flux0',flux0 print 'max_speed', max_speed #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) flux, max_speed = flux_function(-1.0, qr, ql, zl, zr) print flux0 print flux assert allclose(flux + flux0, 0.) #print flux0 #print max_speed #assert allclose(flux0, [-15.3598804, 253.71111111, 0.]) assert allclose(flux0, [-15.3598804, 253.71111111]) 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=2, edge=0) ql = domain.get_conserved_quantities(vol_id=1, vertex=0) qr = domain.get_conserved_quantities(vol_id=0, vertex=1) #flux1, max_speed = flux_function(normal, ql, qr, zl, zr) flux1, max_speed = flux_function(-1.0, ql, qr, zl, zr) #assert allclose(flux1, [2.4098563, 0., 123.04444444]) print 'flux1', flux1 assert allclose(flux1, [2.4098563, 0.]) 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) #ql = domain.get_conserved_quantities(vol_id=1, vertex=2) #qr = domain.get_conserved_quantities(vol_id=0, vertex=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] total_flux = -(flux0+flux1)/domain.areas[1] #assert allclose(total_flux, [-0.68218178, -166.6, -35.93333333]) assert allclose(total_flux, [-0.68218178, -166.6]) domain.compute_fluxes() #assert allclose(total_flux, domain.explicit_update[1,:]) for i, name in enumerate(['stage', '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['stage'].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_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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) 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('stage', [[val0, val0-1, val0-2], # [val1, val1+1, val1], # [val2, val2-2, val2], # [val3-0.5, val3, val3]]) domain.set_quantity('elevation', zl*ones( (4,2) )) domain.set_quantity('stage', [[val0, val0-1], [val1, val1+1], [val2, val2-2], [val3-0.5, val3]]) #Should fail try: domain.check_integrity() except: pass def test_initial_condition(self): """Test that initial condition is output at time == 0 """ from config import g import copy #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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) #Set up for a gradient of (3,0) at mid triangle #def slope(x, y): # return 3*x def slope(x): return 3*x h = 0.1 #def stage(x,y): # return slope(x,y)+h def stage(x): return slope(x)+h domain.set_quantity('elevation', slope) domain.set_quantity('stage', stage) initial_stage = copy.copy(domain.quantities['stage'].vertex_values) domain.set_boundary({'exterior': Reflective_boundary(domain)}) #Evolution for t in domain.evolve(yieldstep = 1.0, finaltime = 2.0): stage = domain.quantities['stage'].vertex_values if t == 0.0: assert allclose(stage, initial_stage) else: assert not allclose(stage, initial_stage) os.remove(domain.filename + '.sww') ##################################################### 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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) def slope(x): return 3*x h = 0.1 #def stage(x,y): # return slope(x,y)+h def stage(x): return slope(x)+h domain.set_quantity('elevation', slope) domain.set_quantity('stage', stage) 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['stage'].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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) #Set up for a gradient of (3,0) at mid triangle def slope(x): return 3*x h = 0.1 #def stage(x,y): # return slope(x,y)+h def stage(x): return slope(x)+h eta = 0.07 domain.set_quantity('elevation', slope) domain.set_quantity('stage', stage) 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['stage'].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['stage'].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['stage'].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['stage'].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('xmomentum', 5) #domain.set_quantity('ymomentum', 4) S = -g * eta**2 * 5 / h**(7.0/3) #S = -g * eta**2 / h**(7.0/3) domain.compute_forcing_terms() assert allclose(domain.quantities['stage'].semi_implicit_update, 0) print '5*S', 5*S #CHECK THIS TEST print domain.quantities['xmomentum'].semi_implicit_update #assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 3*S) assert allclose(domain.quantities['xmomentum'].semi_implicit_update, 5*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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) #Flat surface with 1m of water domain.set_quantity('elevation', 0) domain.set_quantity('stage', 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) S = const * u assert allclose(domain.quantities['stage'].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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) #Flat surface with 1m of water domain.set_quantity('elevation', 0) domain.set_quantity('stage', 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] x = xc #s_vec = speed(t,x,y) #phi_vec = angle(t,x,y) s_vec = speed(t,x) phi_vec = angle(t,x) 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) S = const*u assert allclose(domain.quantities['stage'].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_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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) 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 stage) #domain.set_quantity('elevation', zl*ones( (4,3) )) domain.set_quantity('elevation', zl*ones( (4,2) )) #domain.set_quantity('stage', [[val0, val0-1, val0-2], # [val1, val1+1, val1], # [val2, val2-2, val2], # [val3-0.5, val3, val3]]) domain.set_quantity('stage', [[val0, val0-1], [val1, val1+1], [val2, val2-2], [val3-0.5, val3]]) domain.order = 1 domain.distribute_to_vertices_and_edges() #Check that centroid values were distributed to vertices C = domain.quantities['stage'].centroid_values #for i in range(3): for i in range(2): assert allclose( domain.quantities['stage'].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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) 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('stage', [[val0, val0, val0], # [val1, val1, val1], # [val2, val2, val2], # [val3, val3, val3]]) domain.set_quantity('elevation', [[0,0], [6,0], [6,6], [6,6]]) domain.set_quantity('stage', [[val0, val0], [val1, val1], [val2, val2], [val3, val3]]) E = domain.quantities['elevation'].vertex_values L = domain.quantities['stage'].vertex_values #Check that some stages 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 stages 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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) val0 = 2. val1 = 4. val2 = 8. val3 = 2. domain.set_quantity('stage', [val0, val1, val2, val3], location='centroids') L = domain.quantities['stage'].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] a=0.0 b=2.0 c=4.0 d=6.0 e=8.0 #f=10.0 #points = [a, b, c, d, e, f] points = [a, b, c, d, e] #bac, bce, ecf, dbe #vertices = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4]] #domain = Domain(points, vertices) domain = Domain(points) L = domain.quantities['stage'].vertex_values #def stage(x,y): # return x**2 def stage(x): return x**2 domain.set_quantity('stage', stage, location='centroids') a = domain.quantities['stage'].compute_gradients() print 'a1', a ###assert allclose(a[1], 3.33333334) #assert allclose(b[1], 0.0) domain.order = 1 domain.distribute_to_vertices_and_edges() print 'l',L #assert allclose(L[1], 1.77777778) assert allclose(L[1], 9) domain.order = 2 domain.distribute_to_vertices_and_edges() assert allclose(L[1], [0.57777777, 2.37777778, 2.37777778]) #------------------------------------------------------------- if __name__ == "__main__": suite = unittest.makeSuite(TestCase,'test') runner = unittest.TextTestRunner(verbosity=2) runner.run(suite)