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

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

Rip currents from Rowan Romeyn

File size: 6.0 KB
Line 
1"""Simple water flow example using ANUGA
2"""
3
4#------------------------------------------------------------------------------
5# Import necessary modules
6#------------------------------------------------------------------------------
7from anuga.abstract_2d_finite_volumes.mesh_factory import *
8from anuga.shallow_water import *
9from anuga.shallow_water.shallow_water_domain import *
10from math import *
11from numpy import *
12import numpy
13from matplotlib import *
14from pylab import *
15import csv
16
17#Parameters
18
19filename = "WORKING-RIP-LAB_Expt-Geometry_Triangular_Mesh"
20location_of_shore = 140    #the position along the y axis of the shorefront
21sandbar = 1.2 #height of sandbar
22sealevel = 0 #height of coast above sea level
23steepness = 8000 #period of sandbar - larger number gives smoother slope - longer period
24halfchannelwidth = 5
25bank_slope = 0.1
26simulation_length = 1
27timestep = 1
28
29#------------------------------------------------------------------------------
30# Setup computational domain
31#------------------------------------------------------------------------------
32length = 120
33width = 170
34seafloor_resolution = 60.0 # Resolution: Max area of triangles in the mesh for seafloor
35feature_resolution = 1.0
36beach_resolution = 10.0
37
38sea_boundary_polygon = [[0,0],[length,0],[length,width],[0,width]]
39feature_boundary_polygon = [[0,100],[length,100],[length,150],[0,150]]
40beach_interior_polygon = [[0,150],[length,150],[length,width],[0,width]]
41
42meshname = str(filename)+'.msh'
43
44#interior regions
45feature_regions = [[feature_boundary_polygon, feature_resolution],
46                   [beach_interior_polygon, beach_resolution]]
47
48create_mesh_from_regions(sea_boundary_polygon,
49                         boundary_tags={'bottom': [0],
50                                        'right' : [1],
51                                        'top' :   [2],
52                                        'left':   [3]},
53                         maximum_triangle_area = seafloor_resolution,
54                         filename = meshname,
55                         interior_regions = feature_regions,
56                         use_cache = False,
57                         verbose = False)
58domain = Domain(meshname, use_cache=True, verbose=True)
59domain.set_name(filename)                  # Output name
60print domain.statistics()
61
62#------------------------------------------------------------------------------
63# Setup initial conditions
64#------------------------------------------------------------------------------
65
66def topography(x,y):
67    """Complex topography defined by a function of vectors x and y."""
68
69    #general slope and buildings
70    z=0.05*(y-(location_of_shore))
71   
72    N = len(x)
73    for i in range(N):
74        if y[i] < 25:
75            z[i] = (0.2*(y[i]-25)) + 0.05*(y[i]-(location_of_shore))
76    for i in range(N):
77        if y[i]>150:
78            z[i] = (0.1*(y[i]-150)) + 0.05*(y[i] - (location_of_shore))
79
80    return z
81
82
83def topography3(x,y):
84    z=0*x
85   
86    N = len(x)
87    for i in range(N):
88        if -1*(bank_slope)*x[i] + 112 < y[i] < -1*(bank_slope)*x[i] + 124 and 0<x[i]<((length/2)-halfchannelwidth):
89            z[i] += sandbar*((cos((y[i]-118)/steepness)))
90    for i in range(N):
91        if -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+112) < y[i] < -1*(bank_slope)*(x[i]-(length/2)) + (-1*(bank_slope)*(length/2)+124) and ((length/2)+halfchannelwidth)<x[i]<183:
92            z[i] += sandbar*(cos((y[i]-118)/steepness))
93    return z
94
95domain.set_quantity('elevation', topography)           # elevation is a function
96domain.set_quantity('friction', 0.01)                  # Constant friction
97
98domain.add_quantity('elevation', topography3)
99
100domain.set_quantity('stage', 0)   # Dry initial condition
101
102#------------------------------------------------------------------------------
103# Setup boundary conditions
104#------------------------------------------------------------------------------
105Bi = Dirichlet_boundary([0.4, 0, 0])          # Inflow
106Br = Reflective_boundary(domain)              # Solid reflective wall
107Bo = Dirichlet_boundary([-5, 0, 0])           # Outflow
108
109
110amplitude = 0.4 #amplitude of wave (wave height m)
111period =  5     #wave period (sec)
112
113def wave(t):
114
115    A = amplitude # Amplitude [m] (Wave height)
116    T = period  # Wave period [s]
117
118    if t < 30000000000:
119        return [A*sin(2*pi*t/T) + 1, 0, 0] 
120    else:
121        return [0.0, 0, 0]
122
123Bt = Time_boundary(domain, f=wave)
124
125
126domain.set_boundary({'left': Br, 'right': Br, 'top': Bo, 'bottom': Bt})
127
128#------------------------------------------------------------------------------
129# Evolve system through time
130#------------------------------------------------------------------------------
131
132list1 = [x for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
133
134print list1
135
136u = numpy.zeros(len(list1), 'd')
137v = numpy.zeros(len(list1), 'd')
138
139
140for t in domain.evolve(yieldstep = (timestep), finaltime = (simulation_length)):
141    print domain.timestepping_statistics()
142    S = domain.get_quantity('stage').get_values(interpolation_points=list1)
143    E = domain.get_quantity('elevation').get_values(interpolation_points=list1)
144    depth = S-E
145
146    uh = domain.get_quantity('xmomentum').get_values(interpolation_points=list1)
147    vh = domain.get_quantity('ymomentum').get_values(interpolation_points=list1)
148    u += uh/depth
149    v += vh/depth
150
151n_time_intervals = (simulation_length)/(timestep)
152u_average = u / (n_time_intervals)
153v_average = v / (n_time_intervals)
154
155print "there were", n_time_intervals, "time steps"
156
157print "sum y velocity", v
158print "average y velocity", v_average
159print "sum x velocity", u
160print "average x velocity", u_average
161
162x_output = file("x_velocity.txt", 'w')
163y_output = file("y_velocity.txt", 'w')
164
165print >> x_output, " "
166print >> y_output, " "
167
168print >> x_output, u_average
169print >> y_output, v_average
170
171X = [x[0] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
172Y = [x[1] for x in csv.reader(open('New_gauges.csv','r'),dialect='excel',delimiter=",") ]
173U = u_average.tolist()
174V = v_average.tolist()
175
176print "U = ", U
177print "U has type", type(U)
178
179
180from matplotlib import *
181from pylab import *
182
183figure()
184quiver(X,Y,U,V)
185show()
186
187
Note: See TracBrowser for help on using the repository browser.