source: anuga_work/development/classroom/ripcurrent.py @ 7589

Last change on this file since 7589 was 7589, checked in by ole, 14 years ago

Removed interpolation and brought run time down from 20+ seconds to < 1 second

File size: 8.5 KB
RevLine 
[7581]1"""ANUGA simulation of simple rip current.
2
[7587]3Source: Geometry and wave properties loosely based on those presented in
4OBSERVATIONS OF LABORATORY RIP CURRENTS by Brian K. Sapp,
5School of Civil and Environmental Engineering
6Georgia Institute of Technology
7May 2006
8
[7588]9I will need to make a version which has the exact same geometry as the
10Georgia Tech wavetank if we wish to use a comparison to the results of
11this study as ANUGA validation as i played with the geometry somewhat
12as i completed this model.
[7578]13"""
14
15#------------------------------------------------------------------------------
16# Import necessary modules
17#------------------------------------------------------------------------------
[7579]18from anuga.interface import create_domain_from_regions
19from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
20from anuga.shallow_water.shallow_water_domain import Reflective_boundary
21from anuga.shallow_water.shallow_water_domain import Time_boundary
[7582]22from anuga.shallow_water.data_manager import get_mesh_and_quantities_from_file
[7588]23from pylab import figure, plot, axis, quiver, quiverkey, show, title, axhline
24from pylab import cos, sin, pi
[7578]25import numpy
26import csv
[7579]27import time
[7578]28
29
[7579]30#------------------------------------------------------------------------------
31# Parameters
32#------------------------------------------------------------------------------
[7581]33
[7588]34filename = 'WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh'
[7581]35
[7579]36location_of_shore = 140 # The position along the y axis of the shorefront
37sandbar = 1.2           # Height of sandbar
38sealevel = 0            # Height of coast above sea level
39steepness = 8000        # Period of sandbar -
40                        # larger number gives smoother slope - longer period
[7578]41halfchannelwidth = 5
42bank_slope = 0.1
43simulation_length = 1
44timestep = 1
45
[7579]46
[7578]47#------------------------------------------------------------------------------
48# Setup computational domain
49#------------------------------------------------------------------------------
50length = 120
51width = 170
[7579]52seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh
[7578]53feature_resolution = 1.0
54beach_resolution = 10.0
55
56sea_boundary_polygon = [[0,0],[length,0],[length,width],[0,width]]
57feature_boundary_polygon = [[0,100],[length,100],[length,150],[0,150]]
58beach_interior_polygon = [[0,150],[length,150],[length,width],[0,width]]
59
60meshname = str(filename)+'.msh'
61
[7579]62# Interior regions
[7578]63feature_regions = [[feature_boundary_polygon, feature_resolution],
64                   [beach_interior_polygon, beach_resolution]]
65
[7579]66domain = create_domain_from_regions(sea_boundary_polygon,
67                                    boundary_tags={'bottom': [0],
68                                                   'right' : [1],
69                                                   'top' :   [2],
70                                                   'left':   [3]},
[7583]71                                    maximum_triangle_area=seafloor_resolution,
72                                    mesh_filename=meshname,
73                                    interior_regions=feature_regions,
74                                    use_cache=True,
75                                    verbose=True)
[7579]76                                   
[7581]77domain.set_name(filename) # Output name
[7578]78print domain.statistics()
79
[7579]80
[7578]81#------------------------------------------------------------------------------
82# Setup initial conditions
83#------------------------------------------------------------------------------
84def topography(x,y):
85    """Complex topography defined by a function of vectors x and y."""
86
[7587]87    # General slope, sets the shore at the location defined previously
[7581]88    z=0.05*(y-location_of_shore)
[7587]89
[7588]90    # Steeper slope close to the seaward boundary giving a region of deep water
[7578]91    N = len(x)
92    for i in range(N):
93        if y[i] < 25:
[7581]94            z[i] = 0.2*(y[i]-25) + 0.05*(y[i]-location_of_shore)
[7588]95           
96    # Steeper slope close to the landward boundary, simulating a beach etc
97    # This helps to prevent too much reflection of wave energy off the
98    # landward boundary         
[7578]99    for i in range(N):
100        if y[i]>150:
[7581]101            z[i] = 0.1*(y[i]-150) + 0.05*(y[i]-location_of_shore)
[7578]102
103    return z
104
105
106def topography3(x,y):
107    z=0*x
108   
109    N = len(x)
[7579]110   
[7588]111    # Set up the left hand side of the sandbank
112    # amount which it deviates from parallel with the beach is controlled
113    # by 'bank_slope'
114    # width of the channel (the gap between the two segments of the sandbank)
115    # is controlled by 'halfchannelwidth'
116    # The height of the sandbar is controlled by 'sandbar'
117    # 'steepness' provides control over the slope of the soundbar
118    # (smaller values give a more rounded shape, if too small will produce
119    # peaks and troughs)
[7578]120    for i in range(N):
[7581]121        ymin = -bank_slope*x[i] + 112 
122        ymax = -bank_slope*x[i] + 124
[7579]123        xmin = 0
[7581]124        xmax = length/2-halfchannelwidth
[7579]125        if ymin < y[i] < ymax and xmin < x[i]< xmax:
[7581]126            z[i] += sandbar*cos((y[i]-118)/steepness)
[7579]127   
[7588]128    # Set up the right hand side of the sandbank
129    # changing the sign in y min and y max allows the two halves of the
130    # sandbank to form a v shape
[7578]131    for i in range(N):
[7581]132        ymin = -bank_slope*(x[i]-length/2) - bank_slope*length/2 + 112
133        ymax = -bank_slope*(x[i]-length/2) - bank_slope*length/2 + 124
134        xmin = length/2+halfchannelwidth
[7579]135        xmax = 183
136        if ymin < y[i] < ymax and xmin < x[i] < xmax:
[7581]137            z[i] += sandbar*cos((y[i]-118)/steepness)
[7579]138           
[7578]139    return z
140
[7579]141domain.set_quantity('elevation', topography)  # Apply base elevation function
142domain.add_quantity('elevation', topography3) # Add elevation modification
143domain.set_quantity('friction', 0.01)         # Constant friction
[7588]144domain.set_quantity('stage', 0)               # Constant initial condition at
145                                              # mean sea level
[7578]146
147
148#------------------------------------------------------------------------------
149# Setup boundary conditions
150#------------------------------------------------------------------------------
151Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
152Br = Reflective_boundary(domain)              # Solid reflective wall
153Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
154
155def wave(t):
[7579]156    """Define wave driving the system
157    """
158   
159    A = 0.4  # Amplitude of wave [m] (wave height)
160    T = 5    # Wave period [s]
[7578]161
162    if t < 30000000000:
163        return [A*sin(2*pi*t/T) + 1, 0, 0] 
164    else:
165        return [0.0, 0, 0]
166
167Bt = Time_boundary(domain, f=wave)
168
169
170domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt})
171
[7579]172
[7578]173#------------------------------------------------------------------------------
174# Evolve system through time
175#------------------------------------------------------------------------------
176
[7579]177# Allocate space for velocity values
[7589]178u = numpy.zeros(len(domain))
179v = numpy.zeros(len(domain))
[7578]180
[7579]181t0 = time.time()
182for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length):
[7578]183    print domain.timestepping_statistics()
[7589]184   
185    S = domain.get_quantity('stage').get_values(location='centroids')
186    E = domain.get_quantity('elevation').get_values(location='centroids')
[7578]187    depth = S-E
188
[7589]189    uh = domain.get_quantity('xmomentum').get_values(location='centroids')
190    vh = domain.get_quantity('ymomentum').get_values(location='centroids')
[7578]191    u += uh/depth
192    v += vh/depth
193
[7579]194   
195#------------------------------------------------------------------------------
196# Post processing
[7582]197#------------------------------------------------------------------------------
[7578]198
[7579]199n_time_intervals = simulation_length/timestep
[7581]200print 'There were %i time steps' % n_time_intervals
[7579]201
[7589]202nodes = domain.get_nodes()
203
204X = nodes[:,0]
205Y = nodes[:,1]
[7583]206U = u/n_time_intervals
207V = v/n_time_intervals
[7578]208
[7581]209print 'Computation took %.2f seconds' % (time.time()-t0)
[7578]210
[7588]211key_auto_length = (max(V))/5 # Make the key vector a sensible length not
212                             # sure how to label it with the correct value
213                             # though
[7587]214
215
[7578]216figure()
[7587]217Q = quiver(X,Y,U,V)
[7588]218qk = quiverkey(Q, 0.8, 0.05, key_auto_length, r'$unknown \frac{m}{s}$',
219               labelpos='E', # Need to get the label to show the value of key_auto_length
220               coordinates='figure', 
221               fontproperties={'weight': 'bold'})
222               
[7587]223axis([-10,length + 10, -10, width +10])
224title('Simulation of a Rip-Current, Average Velocity Vector Field')
225
226axhline(y=25,color='b')
227axhline(y=(location_of_shore),color='r')
228
[7588]229x1 = numpy.arange(0,55,1)
[7587]230y1 = -(bank_slope)*x1 + 112
231y12 = -(bank_slope)*x1 + 124
232
[7588]233x2 = numpy.arange(65,length,1)
[7587]234y2 = -(bank_slope)*x2 + 112
235y22 = -(bank_slope)*x2 + 124
236
237plot(x1,y1,x1,y12,x2,y2,x2,y22,color='g')
[7578]238show()
239
240
Note: See TracBrowser for help on using the repository browser.