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

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

Cleanup of school exercise. Optimisation to follow

File size: 8.3 KB
RevLine 
[7578]1"""Simple water flow example using ANUGA
2"""
3
4#------------------------------------------------------------------------------
5# Import necessary modules
6#------------------------------------------------------------------------------
[7579]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
[7578]12import numpy
13import csv
[7579]14import time
[7578]15
16
[7579]17#------------------------------------------------------------------------------
18# Parameters
19#------------------------------------------------------------------------------
[7578]20filename = "WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh"
[7579]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
[7578]26halfchannelwidth = 5
27bank_slope = 0.1
28simulation_length = 1
29timestep = 1
30
[7579]31
[7578]32#------------------------------------------------------------------------------
33# Setup computational domain
34#------------------------------------------------------------------------------
35length = 120
36width = 170
[7579]37seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh
[7578]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
[7579]47# Interior regions
[7578]48feature_regions = [[feature_boundary_polygon, feature_resolution],
49                   [beach_interior_polygon, beach_resolution]]
50
[7579]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                                   
[7578]62domain.set_name(filename)                  # Output name
63print domain.statistics()
64
[7579]65
[7578]66#------------------------------------------------------------------------------
67# Setup initial conditions
68#------------------------------------------------------------------------------
69def topography(x,y):
70    """Complex topography defined by a function of vectors x and y."""
71
[7579]72    # General slope and buildings
[7578]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:
[7579]81            z[i] = (0.1*(y[i]-150)) + 0.05*(y[i]-(location_of_shore))
[7578]82
83    return z
84
85
86def topography3(x,y):
87    z=0*x
88   
89    N = len(x)
[7579]90   
91    # It would be great with a comment about what this does
[7578]92    for i in range(N):
[7579]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:
[7578]98            z[i] += sandbar*((cos((y[i]-118)/steepness)))
[7579]99   
100    # It would be great with a comment about what this does       
[7578]101    for i in range(N):
[7579]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:
[7578]107            z[i] += sandbar*(cos((y[i]-118)/steepness))
[7579]108           
[7578]109    return z
110
[7579]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
[7578]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):
[7579]125    """Define wave driving the system
126    """
127   
128    A = 0.4  # Amplitude of wave [m] (wave height)
129    T = 5    # Wave period [s]
[7578]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
[7579]141
[7578]142#------------------------------------------------------------------------------
143# Evolve system through time
144#------------------------------------------------------------------------------
145
[7579]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                               
[7578]154
[7579]155# Allocate space for velocity values
156u = numpy.zeros(number_of_gauges)
157v = numpy.zeros(number_of_gauges)
[7578]158
[7579]159t0 = time.time()
160for t in domain.evolve(yieldstep = timestep, finaltime = simulation_length):
[7578]161    print domain.timestepping_statistics()
[7579]162    S = domain.get_quantity('stage').get_values(interpolation_points=gauges)
163    E = domain.get_quantity('elevation').get_values(interpolation_points=gauges)
[7578]164    depth = S-E
165
[7579]166    uh = domain.get_quantity('xmomentum').get_values(interpolation_points=gauges)
167    vh = domain.get_quantity('ymomentum').get_values(interpolation_points=gauges)
[7578]168    u += uh/depth
169    v += vh/depth
[7579]170print 'Evolution took %.2f seconds' % (time.time()-t0)
[7578]171
[7579]172   
173#------------------------------------------------------------------------------
174# Post processing
175#------------------------------------------------------------------------------   
[7578]176
[7579]177# Use this
[7578]178
[7579]179#     from anuga.fit_interpolate.interpolate import Interpolation_function
[7578]180
[7579]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
[7578]186
[7579]187#     # Find all intersections and associated triangles.
188#     segments = mesh.get_intersecting_segments(polyline, verbose=verbose)
[7578]189
[7579]190#     # Get midpoints
191#     interpolation_points = segment_midpoints(segments)
[7578]192
[7579]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                             
[7578]237U = u_average.tolist()
238V = v_average.tolist()
239
[7579]240#print "U = ", U
241#print "U has type", type(U)
[7578]242
243
244figure()
245quiver(X,Y,U,V)
246show()
247
248
Note: See TracBrowser for help on using the repository browser.