"""ANUGA simulation of simple rip current. Source: Geometry and wave properties loosely based on those presented in OBSERVATIONS OF LABORATORY RIP CURRENTS by Brian K. Sapp, School of Civil and Environmental Engineering Georgia Institute of Technology May 2006 I will need to make a version which has the exact same geometry as the Georgia Tech wavetank if we wish to use a comparison to the results of this study as ANUGA validation as i played with the geometry somewhat as i completed this model. """ #------------------------------------------------------------------------------ # Import necessary modules #------------------------------------------------------------------------------ import anuga from pylab import figure, plot, axis, quiver, quiverkey, show, title, axhline from pylab import cos, sin, pi import numpy import csv import time #------------------------------------------------------------------------------ # Parameters #------------------------------------------------------------------------------ filename = 'WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh' location_of_shore = 140 # The position along the y axis of the shorefront sandbar = 1.2 # Height of sandbar sealevel = 0 # Height of coast above sea level steepness = 8000 # Period of sandbar - # larger number gives smoother slope - longer period halfchannelwidth = 5 bank_slope = 0.1 simulation_length = 60 timestep = 1 #------------------------------------------------------------------------------ # Setup computational domain #------------------------------------------------------------------------------ length = 120 width = 170 seafloor_resolution = 20.0 # Resolution: Max area of triangles in the mesh feature_resolution = 1.0 beach_resolution = 10.0 sea_boundary_polygon = [[0,0],[length,0],[length,width],[0,width]] feature_boundary_polygon = [[19,99],[length/2+1,99],[length/2+1,151],[0,151]] hole_boundary_polygon = [[20,100],[length/2,100],[length/2,150],[20,150]] beach_interior_polygon = [[0,150],[length,150],[length,width],[0,width]] meshname = str(filename)+'.msh' # Interior regions feature_regions = [[feature_boundary_polygon, feature_resolution], [beach_interior_polygon, beach_resolution]] domain = anuga.create_domain_from_regions(sea_boundary_polygon, boundary_tags={'bottom': [0], 'right' : [1], 'top' : [2], 'left': [3]}, maximum_triangle_area=seafloor_resolution, mesh_filename=meshname, interior_regions=feature_regions, interior_holes=[hole_boundary_polygon], use_cache=True, verbose=True) domain.set_name(filename) # Output name print domain.statistics() #------------------------------------------------------------------------------ # Setup initial conditions #------------------------------------------------------------------------------ def topography(x,y): """Complex topography defined by a function of vectors x and y.""" # General slope, sets the shore at the location defined previously z=0.05*(y-location_of_shore) # Steeper slope close to the seaward boundary giving a region of deep water N = len(x) for i in range(N): if y[i] < 25: z[i] = 0.2*(y[i]-25) + 0.05*(y[i]-location_of_shore) # Steeper slope close to the landward boundary, simulating a beach etc # This helps to prevent too much reflection of wave energy off the # landward boundary for i in range(N): if y[i]>150: z[i] = 0.1*(y[i]-150) + 0.05*(y[i]-location_of_shore) return z def topography3(x,y): z=0*x N = len(x) # Set up the left hand side of the sandbank # amount which it deviates from parallel with the beach is controlled # by 'bank_slope' # width of the channel (the gap between the two segments of the sandbank) # is controlled by 'halfchannelwidth' # The height of the sandbar is controlled by 'sandbar' # 'steepness' provides control over the slope of the soundbar # (smaller values give a more rounded shape, if too small will produce # peaks and troughs) for i in range(N): ymin = -bank_slope*x[i] + 112 ymax = -bank_slope*x[i] + 124 xmin = 0 xmax = length/2-halfchannelwidth if ymin < y[i] < ymax and xmin < x[i]< xmax: z[i] += sandbar*cos((y[i]-118)/steepness) # Set up the right hand side of the sandbank # changing the sign in y min and y max allows the two halves of the # sandbank to form a v shape for i in range(N): ymin = -bank_slope*(x[i]-length/2) - bank_slope*length/2 + 112 ymax = -bank_slope*(x[i]-length/2) - bank_slope*length/2 + 124 xmin = length/2+halfchannelwidth xmax = 183 if ymin < y[i] < ymax and xmin < x[i] < xmax: z[i] += sandbar*cos((y[i]-118)/steepness) return z domain.set_quantity('elevation', topography) # Apply base elevation function domain.add_quantity('elevation', topography3) # Add elevation modification domain.set_quantity('friction', 0.01) # Constant friction domain.set_quantity('stage', 0) # Constant initial condition at # mean sea level #------------------------------------------------------------------------------ # Setup boundary conditions #------------------------------------------------------------------------------ Bi = anuga.Dirichlet_boundary([0.4, 0, 0]) # Inflow Br = anuga.Reflective_boundary(domain) # Solid reflective wall Bo = anuga.Dirichlet_boundary([-5, 0, 0]) # Outflow def wave(t): """Define wave driving the system """ A = 0.4 # Amplitude of wave [m] (wave height) T = 1 # Wave period [s] if t < 30000000000: return [A*sin(2*pi*t/T) + 1, 0, 0] else: return [0.0, 0, 0] Bt = anuga.Time_boundary(domain, f=wave) domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt, 'exterior': Br}) #------------------------------------------------------------------------------ # Evolve system through time #------------------------------------------------------------------------------ # Allocate space for velocity values u = numpy.zeros(len(domain)) v = numpy.zeros(len(domain)) t0 = time.time() for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length): print domain.timestepping_statistics() S = domain.get_quantity('stage').get_values(location='centroids') E = domain.get_quantity('elevation').get_values(location='centroids') depth = S-E uh = domain.get_quantity('xmomentum').get_values(location='centroids') vh = domain.get_quantity('ymomentum').get_values(location='centroids') u += uh/depth v += vh/depth #------------------------------------------------------------------------------ # Post processing #------------------------------------------------------------------------------ n_time_intervals = simulation_length/timestep print 'There were %i time steps' % n_time_intervals nodes = domain.get_nodes() X = nodes[:,0] Y = nodes[:,1] U = u/n_time_intervals V = v/n_time_intervals print 'Computation took %.2f seconds' % (time.time()-t0) key_auto_length = (max(V))/5 # Make the key vector a sensible length not # sure how to label it with the correct value # though figure() Q = quiver(X,Y,U,V) qk = quiverkey(Q, 0.8, 0.05, key_auto_length, r'$unknown \frac{m}{s}$', labelpos='E', # Need to get the label to show the value of key_auto_length coordinates='figure', fontproperties={'weight': 'bold'}) axis([-10,length + 10, -10, width +10]) title('Simulation of a Rip-Current, Average Velocity Vector Field') axhline(y=25,color='b') axhline(y=(location_of_shore),color='r') x1 = numpy.arange(0,55,1) y1 = -(bank_slope)*x1 + 112 y12 = -(bank_slope)*x1 + 124 x2 = numpy.arange(65,length,1) y2 = -(bank_slope)*x2 + 112 y22 = -(bank_slope)*x2 + 124 plot(x1,y1,x1,y12,x2,y2,x2,y22,color='g') show()