from math import sqrt, pi from shallow_water_1d import * from Numeric import allclose, array, zeros, ones, Float, take from config import g, epsilon def test_sqrt(): for i in range (80000): a = sqrt(4356) b = sqrt(2031) def test_flux1(): #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) print 'flux', flux def test_compute_fluxes0(): #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_1d_solution_I(): print "TEST 1D-SOLUTION I" L = 2000.0 # Length of channel (m) N = 100 # Number of compuational cells cell_len = L/N # Origin = 0.0 points = zeros(N+1,Float) for i in range(N+1): points[i] = i*cell_len domain = Domain(points) domain.order = 2 def stage(x): for i in range(len(x)): if x[i]<=1000.0: x[i] = 10.0 else: x[i] = 5.0 return x domain.set_quantity('stage', stage) #L = domain.quantities['stage'].vertex_values #print "Initial Stage" #print L domain.set_boundary({'exterior': Reflective_boundary(domain)}) import time t0 = time.time() yieldstep = 1.0 finaltime = 50.0 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): #xmom = domain.quantities['xmomentum'] #xmom = xmom.centroid_values #stage = domain.quantities['stage'] #stage = stage.centroid_values #print 'stage', stage #print 'xmom', xmom #for i in range N #u = xmom/stage pass print 'That took %.2f seconds' %(time.time()-t0) #L = domain.quantities['stage'].vertex_values #print "Final Stage" #print L C = domain.quantities['stage'].vertex_values #print C f = file('test_solution_Ix.out', 'w') for i in range(N): f.write(str(C[i,1])) f.write("\n") f.close