source: trunk/anuga_work/development/classroom/ripcurrent.py @ 7877

Last change on this file since 7877 was 7877, checked in by hudson, 14 years ago

Moved all development files into trunk.

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