source: anuga_work/development/classroom/WORKING-RipCurrent.py @ 7579

Last change on this file since 7579 was 7579, checked in by ole, 13 years ago

Cleanup of school exercise. Optimisation to follow

File size: 8.3 KB
Line 
1"""Simple water flow example using ANUGA
2"""
3
4#------------------------------------------------------------------------------
5# Import necessary modules
6#------------------------------------------------------------------------------
7from anuga.interface import create_domain_from_regions
8from anuga.shallow_water.shallow_water_domain import Dirichlet_boundary
9from anuga.shallow_water.shallow_water_domain import Reflective_boundary
10from anuga.shallow_water.shallow_water_domain import Time_boundary
11from pylab import figure, quiver, show, cos, sin, pi
12import numpy
13import csv
14import time
15
16
17#------------------------------------------------------------------------------
18# Parameters
19#------------------------------------------------------------------------------
20filename = "WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh"
21location_of_shore = 140 # The position along the y axis of the shorefront
22sandbar = 1.2           # Height of sandbar
23sealevel = 0            # Height of coast above sea level
24steepness = 8000        # Period of sandbar -
25                        # larger number gives smoother slope - longer period
26halfchannelwidth = 5
27bank_slope = 0.1
28simulation_length = 1
29timestep = 1
30
31
32#------------------------------------------------------------------------------
33# Setup computational domain
34#------------------------------------------------------------------------------
35length = 120
36width = 170
37seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh
38feature_resolution = 1.0
39beach_resolution = 10.0
40
41sea_boundary_polygon = [[0,0],[length,0],[length,width],[0,width]]
42feature_boundary_polygon = [[0,100],[length,100],[length,150],[0,150]]
43beach_interior_polygon = [[0,150],[length,150],[length,width],[0,width]]
44
45meshname = str(filename)+'.msh'
46
47# Interior regions
48feature_regions = [[feature_boundary_polygon, feature_resolution],
49                   [beach_interior_polygon, beach_resolution]]
50
51domain = create_domain_from_regions(sea_boundary_polygon,
52                                    boundary_tags={'bottom': [0],
53                                                   'right' : [1],
54                                                   'top' :   [2],
55                                                   'left':   [3]},
56                                    maximum_triangle_area = seafloor_resolution,
57                                    mesh_filename = meshname,
58                                    interior_regions = feature_regions,
59                                    use_cache = True,
60                                    verbose = True)
61                                   
62domain.set_name(filename)                  # Output name
63print domain.statistics()
64
65
66#------------------------------------------------------------------------------
67# Setup initial conditions
68#------------------------------------------------------------------------------
69def topography(x,y):
70    """Complex topography defined by a function of vectors x and y."""
71
72    # General slope and buildings
73    z=0.05*(y-(location_of_shore))
74   
75    N = len(x)
76    for i in range(N):
77        if y[i] < 25:
78            z[i] = (0.2*(y[i]-25)) + 0.05*(y[i]-(location_of_shore))
79    for i in range(N):
80        if y[i]>150:
81            z[i] = (0.1*(y[i]-150)) + 0.05*(y[i]-(location_of_shore))
82
83    return z
84
85
86def topography3(x,y):
87    z=0*x
88   
89    N = len(x)
90   
91    # It would be great with a comment about what this does
92    for i in range(N):
93        ymin = -1*(bank_slope)*x[i] + 112 
94        ymax = -1*(bank_slope)*x[i] + 124
95        xmin = 0
96        xmax = (length/2)-halfchannelwidth
97        if ymin < y[i] < ymax and xmin < x[i]< xmax:
98            z[i] += sandbar*((cos((y[i]-118)/steepness)))
99   
100    # It would be great with a comment about what this does       
101    for i in range(N):
102        ymin = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112)
103        ymax = -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124)
104        xmin = (length/2)+halfchannelwidth
105        xmax = 183
106        if ymin < y[i] < ymax and xmin < x[i] < xmax:
107            z[i] += sandbar*(cos((y[i]-118)/steepness))
108           
109    return z
110
111domain.set_quantity('elevation', topography)  # Apply base elevation function
112domain.add_quantity('elevation', topography3) # Add elevation modification
113domain.set_quantity('friction', 0.01)         # Constant friction
114domain.set_quantity('stage', 0)               # Constant initial condition at mean sea level
115
116
117#------------------------------------------------------------------------------
118# Setup boundary conditions
119#------------------------------------------------------------------------------
120Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
121Br = Reflective_boundary(domain)              # Solid reflective wall
122Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
123
124def wave(t):
125    """Define wave driving the system
126    """
127   
128    A = 0.4  # Amplitude of wave [m] (wave height)
129    T = 5    # Wave period [s]
130
131    if t < 30000000000:
132        return [A*sin(2*pi*t/T) + 1, 0, 0] 
133    else:
134        return [0.0, 0, 0]
135
136Bt = Time_boundary(domain, f=wave)
137
138
139domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt})
140
141
142#------------------------------------------------------------------------------
143# Evolve system through time
144#------------------------------------------------------------------------------
145
146# Read in gauge locations for interpolation and convert them to floats
147gauge_file = open('New_gauges.csv', 'r')
148G = [(float(x[0]), float(x[1])) for x in csv.reader(gauge_file,
149                                                    dialect='excel', 
150                                                    delimiter=',')]
151gauges = numpy.array(G)
152number_of_gauges = len(gauges)
153print gauges                               
154
155# Allocate space for velocity values
156u = numpy.zeros(number_of_gauges)
157v = numpy.zeros(number_of_gauges)
158
159t0 = time.time()
160for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length):
161    print domain.timestepping_statistics()
162    S = domain.get_quantity('stage').get_values(interpolation_points=gauges)
163    E = domain.get_quantity('elevation').get_values(interpolation_points=gauges)
164    depth = S-E
165
166    uh = domain.get_quantity('xmomentum').get_values(interpolation_points=gauges)
167    vh = domain.get_quantity('ymomentum').get_values(interpolation_points=gauges)
168    u += uh/depth
169    v += vh/depth
170print 'Evolution took %.2f seconds' % (time.time()-t0)
171
172   
173#------------------------------------------------------------------------------
174# Post processing
175#------------------------------------------------------------------------------   
176
177# Use this
178
179#     from anuga.fit_interpolate.interpolate import Interpolation_function
180
181#     # Get mesh and quantities from sww file
182#     X = get_mesh_and_quantities_from_file(filename,
183#                                           quantities=quantity_names,
184#                                           verbose=verbose)
185#     mesh, quantities, time = X
186
187#     # Find all intersections and associated triangles.
188#     segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
189
190#     # Get midpoints
191#     interpolation_points = segment_midpoints(segments)
192
193#     # Interpolate
194#     if verbose:
195#         log.critical('Interpolating - total number of interpolation points = %d'
196#                      % len(interpolation_points))
197
198#     I = Interpolation_function(time,
199#                                quantities,
200#                                quantity_names=quantity_names,
201#                                vertex_coordinates=mesh.nodes,
202#                                triangles=mesh.triangles,
203#                                interpolation_points=interpolation_points,
204#                                verbose=verbose)
205
206#     return segments, I
207
208
209
210
211
212
213n_time_intervals = simulation_length/timestep
214u_average = u/n_time_intervals
215v_average = v/n_time_intervals
216
217#print "there were", n_time_intervals, "time steps"
218
219#print "sum y velocity", v
220#print "average y velocity", v_average
221#print "sum x velocity", u
222#print "average x velocity", u_average
223
224x_output = file('x_velocity.txt', 'w')
225y_output = file('y_velocity.txt', 'w')
226
227#print >> x_output, " "
228#print >> y_output, " "
229
230#print >> x_output, u_average
231#print >> y_output, v_average
232
233
234X = gauges[:,0] 
235Y = gauges[:,1] 
236                             
237U = u_average.tolist()
238V = v_average.tolist()
239
240#print "U = ", U
241#print "U has type", type(U)
242
243
244figure()
245quiver(X,Y,U,V)
246show()
247
248
Note: See TracBrowser for help on using the repository browser.