""" Example of shallow water wave equation consisting of an asymetrical converging channel. Copyright 2004 Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray Geoscience Australia, 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 # #Were these used? #import visualise2_chris as visualise #import Image, ImageGrab import sys from os import sep sys.path.append('..'+sep+'pyvolution') from shallow_water import Domain, Constant_height from shallow_water import Transmissive_boundary, Reflective_boundary,\ Dirichlet_boundary , Transmissive_Momentum_Set_Stage_boundary from math import sqrt, cos, sin, pi from mesh_factory import oblique_cross ###################### # Domain # leny = 30. lenx = 40. f = 5 n = 60*f m = 80*f theta = 25 h_bc = 1.5 p_bc = 5.0 points, elements, boundary = oblique_cross(m, n, lenx, leny, theta = theta) domain = Domain(points, elements, boundary) print 'Number of Elements = ',domain.number_of_elements # Order of solver domain.default_order=2 domain.beta_w = 0.8 # Store output domain.store=True # Provide file name for storing output domain.filename="oblique_cross_%g_%g_%g_%g_%g" % (n,m,theta,h_bc,p_bc) print domain.filename # Output format domain.format="sww" #NET.CDF binary format # "dat" for ASCII # Visualization smoothing domain.smooth=False domain.visualise=False ####################### #Bed-slope and friction class ConstantFunctionT: def __init__(self,value): self.value = value def __call__(self,t): return self.value shock_hb = ConstantFunctionT(h_bc) shock_hh = ConstantFunctionT(h_bc) shock_pb = ConstantFunctionT(p_bc) shock_pp = ConstantFunctionT(p_bc) def x_slope(x, y): return 0*x def shock_h(x,y): n = x.shape[0] w = 0*x for i in range(n): if x[i]<0: w[i] = shock_hh(0.0) else: w[i] = 1.0 return w def shock_p(x,y): n = x.shape[0] w = 0*x for i in range(n): if x[i]<0: w[i] = shock_pp(0.0) else: w[i] = 0.0 return w domain.set_quantity('elevation', x_slope) domain.set_quantity('friction', 0.0) ###################### # Boundary conditions # R = Reflective_boundary(domain) T = Transmissive_boundary(domain) D = Dirichlet_boundary([shock_hb(0.0) , shock_pb(0.0), 0.0]) Bts = Transmissive_Momentum_Set_Stage_boundary(domain, shock_hb) domain.set_boundary({'left': D, 'right': T, 'top': R, 'bottom': R}) ###################### #Initial condition domain.set_quantity('stage', shock_h ) domain.set_quantity('xmomentum',shock_p ) class SetValueWhere: def __init__(self,quantity,value=0,xrange=0.0): self.quantity = quantity self.value = value self.xrange = xrange def __call__(self,x,y): n = x.shape[0] w = self.quantity.centroid_values for i in range(n): if x[i] tclean-0.09 and t < tclean+0.01: # print 'Cleaning up Initial profile' # setstage = SetValueWhere(quantity=Stage,value=vstage[0],xrange=10) # setxmom = SetValueWhere(quantity=Xmom,value=vxmom[0],xrange=10) # Stage.set_values(setstage,location='centroids') # Xmom.set_values(setxmom,location='centroids') # Dnew = Dirichlet_boundary([vstage[0] , vxmom[0], 0.0]) # domain.set_boundary({'left': Dnew, 'right': T, 'top': R, 'bottom': R}) print 'That took %.2f seconds' %(time.time()-t0) #FIXME: Compute average water depth on either side of shock and compare #to expected values. And also Froude numbers. #print "saving file?" #im = ImageGrab.grab() #im.save("ccube.eps")