"""Example of shallow water wave equation analytical solution consists of a flat water surface profile in a parabolic basin with linear friction. The analytical solution was derived by Sampson in 2002. Copyright 2004 Christopher Zoppou, Stephen Roberts ANU Specific methods pertaining to the 2D shallow water equation are imported from shallow_water for use with the generic finite volume framework Conserved quantities are h, uh and vh stored as elements 0, 1 and 2 in the numerical vector named conserved_quantities. """ ###################### # Module imports # import sys from os import sep sys.path.append('..'+sep+'pyvolution') from anuga.pyvolution.shallow_water import Domain, Transmissive_boundary, Reflective_boundary,\ Dirichlet_boundary, gravity, linear_friction from math import sqrt, cos, sin, pi, exp from anuga.pyvolution.mesh_factory import rectangular_cross from anuga.pyvolution.quantity import Quantity from anuga.utilities.polygon import inside_polygon from Numeric import asarray from least_squares import Interpolation #------------------------------- # Domain n = 100 m = 100 lenx = 10000.0 leny = 10000.0 origin = (-5000.0, -5000.0) points, elements, boundary = rectangular_cross(m, n, lenx, leny, origin) domain = Domain(points, elements, boundary) #---------------- # Order of scheme domain.default_order = 1 domain.smooth = True domain.quantities['linear_friction'] = Quantity(domain) #--------------------------------------------------------------------------- # Reconstruct forcing terms with linear friction instead of manning friction domain.forcing_terms = [] domain.forcing_terms.append(gravity) domain.forcing_terms.append(linear_friction) print domain.forcing_terms #------------------------------------- # Provide file name for storing output domain.store = True domain.format = 'sww' domain.set_name('sampson_first_order_cross') print 'Number of triangles = ', len(domain) #----------------------------------------- #Reduction operation for get_vertex_values from anuga.utilities.numerical_tools import mean domain.reduction = mean #domain.reduction = min #Looks better near steep slopes #----------------- #Initial condition # print 'Initial condition' t = 0.0 h0 = 10. a = 3000. g = 9.81 tau =0.001 B = 5 A = 0 p = sqrt(8*g*h0/a/a) s = sqrt(p*p-tau*tau)/2 t = 0. #------------------------------ #Set bed-elevation and friction def x_slope(x,y): n = x.shape[0] z = 0*x for i in range(n): z[i] = h0 - h0*(1.0 -x[i]*x[i]/a/a - y[i]*y[i]/a/a) return z domain.set_quantity('elevation', x_slope) domain.set_quantity('linear_friction', tau) #------------------- #Set the water stage def stage(x,y): z = x_slope(x,y) n = x.shape[0] h = 0*x for i in range(n): h[i] = h0-B*B*exp(-tau*t)/2/g-1/g*(exp(-tau*t/2)*(B*s*cos(s*t) \ +tau*B/2*sin(s*t)))*x[i] \ -1/g*(exp(-tau*t/2)*(B*s*sin(s*t)-tau*B/2*cos(s*t)))*y[i] if h[i] < z[i]: h[i] = z[i] return h domain.set_quantity('stage', stage) #--------- # Boundary print 'Boundary conditions' R = Reflective_boundary(domain) T = Transmissive_boundary(domain) D = Dirichlet_boundary([0.0, 0.0, 0.0]) domain.set_boundary({'left': D, 'right': D, 'top': D, 'bottom': D}) #--------------------------------------------- # Find triangle that contains the point points # and print to file points = [0.,0.] for n in range(len(domain.triangles)): t = domain.triangles[n] tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] if inside_polygon(points,tri): print 'Point is within triangle with vertices '+'%s'%tri n_point = n print 'n_point = ',n_point t = domain.triangles[n_point] tri = [domain.coordinates[t[0]],domain.coordinates[t[1]],domain.coordinates[t[2]]] filename=domain.get_name() file = open(filename,'w') #---------- # Evolution import time t0 = time.time() Stage = domain.quantities['stage'] Xmomentum = domain.quantities['xmomentum'] Ymomentum = domain.quantities['ymomentum'] for t in domain.evolve(yieldstep = 20.0, finaltime = 10000.0 ): domain.write_time() tri_array = asarray(tri) t_array = asarray([[0,1,2]]) interp = Interpolation(tri_array,t_array,[points]) stage = Stage.get_values(location='centroids',indices=[n_point])[0] xmomentum = Xmomentum.get_values(location='centroids',indices=[n_point])[0] ymomentum = Ymomentum.get_values(location='centroids',indices=[n_point])[0] file.write( '%10.6f %10.6f %10.6f %10.6f\n'%(t,stage,xmomentum,ymomentum) ) file.close() print 'That took %.2f seconds' %(time.time()-t0)