source: anuga_core/source/anuga/examples/runup_convergence.py @ 3618

Last change on this file since 3618 was 3618, checked in by ole, 18 years ago

Working meeting with Rosh

File size: 8.6 KB
Line 
1"""Simple water flow example using ANUGA
2
3Water driven up a linear slope and time varying boundary,
4similar to a beach environment.
5
6The study area is discretised as a regular triangular grid 100m x 100m
7"""
8
9
10#------------------------------------------------------------------------------
11# Import necessary modules
12#------------------------------------------------------------------------------
13
14import sys
15
16from anuga.pmesh.mesh_interface import create_mesh_from_regions
17from abstract_2d_finite_volumes.mesh_factory import rectangular_cross
18from anuga.config import g
19from anuga.shallow_water import Domain
20from anuga.shallow_water import Reflective_boundary
21from anuga.shallow_water import Dirichlet_boundary
22from anuga.shallow_water import Time_boundary
23from anuga.shallow_water import Transmissive_Momentum_Set_Stage_boundary
24from abstract_2d_finite_volumes.util import file_function
25from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis, legend, grid, hold
26
27
28
29#------------------------------------------------------------------------------
30# Model constants
31
32slope = -0.02       # 1:50 Slope, reaches h=20m 1000m from western bndry, and h=0 (coast) at 300m
33highest_point = 6   # Highest elevation (m)
34sea_level = 0       # Mean sea level
35min_elevation = -20 # Lowest elevation (elevation of offshore flat part)
36offshore_depth=sea_level-min_elevation # offshore water depth
37amplitude = 0.5       # Solitary wave height H
38normalized_amplitude = amplitude/offshore_depth
39simulation_name = 'runup_convergence'   
40
41
42# Basin dimensions (m)
43west = 0          # left boundary
44east = 1500       # right boundary
45south = 0         # lower boundary
46north = 100       # upper boundary
47
48
49#------------------------------------------------------------------------------
50# Setup computational domain all units in meters
51#------------------------------------------------------------------------------
52
53# Structured mesh
54dx = 20           # Resolution: Length of subdivisions on x axis (length)
55dy = 20           # Resolution: Length of subdivisions on y axis (width)
56
57length = east-west
58width = north-south
59points, vertices, boundary = rectangular_cross(length/dx, width/dy, len1=length, len2=width,
60                                               origin = (west, south)) 
61
62domain = Domain(points, vertices, boundary) # Create domain
63
64
65
66# Unstructured mesh
67polygon = [[east,north],[west,north],[west,south],[east,south]]
68interior_polygon = [[400,north-10],[west+10,north-10],[west+10,south+10],[400,south+10]]
69meshname = simulation_name + '.msh'
70create_mesh_from_regions(polygon,
71                         boundary_tags={'top': [0], 'left': [1], 'bottom': [2], 'right': [3]},
72                         maximum_triangle_area=dx*dy/4, # Triangle area commensurate with structured mesh
73                         filename=meshname,
74                         interior_regions=[[interior_polygon,dx*dy/32]])
75domain = Domain(meshname, use_cache=True, verbose = True)
76
77
78
79domain.set_name(simulation_name)
80
81
82#------------------------------------------------------------------------------
83# Setup initial conditions
84#------------------------------------------------------------------------------
85
86#def topography(x,y):
87#    return slope*x+highest_point  # Return linear bed slope bathymetry as vector
88
89def topography(x,y):
90    """Two part topography - slope and flat part
91    """
92
93    from Numeric import zeros, Float
94
95    z = zeros(len(x), Float) # Allocate space for return vector
96    for i in range(len(x)):
97
98        z[i] = slope*x[i]+highest_point # Linear bed slope bathymetry
99
100        if z[i] < min_elevation:        # Limit depth
101            z[i] = min_elevation
102
103    return z       
104       
105
106       
107
108domain.set_quantity('elevation', topography) # Use function for elevation
109domain.set_quantity('friction', 0.0 )        # Constant friction
110domain.set_quantity('stage', sea_level)      # Constant initial stage
111
112
113#------------------------------------------------------------------------------
114# Setup boundary conditions
115#------------------------------------------------------------------------------
116
117from math import sin, pi, cosh, sqrt
118Br = Reflective_boundary(domain)      # Solid reflective wall
119Bd = Dirichlet_boundary([0.,0.,0.])   # Constant boundary values
120
121
122def waveform(t): 
123    return sea_level + amplitude/cosh(((t-50)/offshore_depth)*(0.75*g*amplitude)**0.5)**2
124
125Bw = Time_boundary(domain=domain,     # Time dependent boundary 
126                   f=lambda t: [waveform(t), 0.0, 0.0])
127
128# Time dependent boundary for stage, where momentum is set automatically
129Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
130
131
132# Associate boundary tags with boundary objects
133domain.set_boundary({'left': Br, 'right': Bts, 'top': Br, 'bottom': Br})
134
135
136#------------------------------------------------------------------------------
137# Evolve system through time
138#------------------------------------------------------------------------------
139
140stagestep = []
141for t in domain.evolve(yieldstep = 1, finaltime = 300):
142    domain.write_time()
143
144    # Let's find the maximum runup here working directly with the quantities,
145    # and stop when it has been detected.
146
147    # 1 Find coastline as x where z==0
148    # 2 Workout travel time to coastline
149    # 3 Find min x where h>0 over all t.
150    # 4 Perhaps do this across a number of ys
151   
152
153#-----------------------------------------------------------------------------
154# Interrogate solution
155#-----------------------------------------------------------------------------
156
157
158# Define line of gauges through center of domain
159def gauge_line(west,east,north,south):
160    from Numeric import arange
161    x_vector = arange(west,600, 10) # Gauges every 1 meter from west to 600m from western bdry
162    y = (north+south)/2.
163
164    gauges = []
165    for x in x_vector:
166        gauges.append([x,y])
167       
168    return gauges, x_vector
169
170gauges, x_vector = gauge_line(west,east,north,south)
171
172# Obtain interpolated timeseries at gauges
173f = file_function(domain.get_name()+'.sww',
174                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
175                  interpolation_points = gauges,
176                  verbose = True,
177                  use_cache = True)
178
179
180# Find runup distance from western boundary through a linear search
181max_stage = []
182min_stage = []
183runup_point = west
184coastline = east       
185for k, g in enumerate(gauges):
186    z = f(0, point_id=k)[1] # Elevation
187
188    min_w = sys.maxint
189    max_w = -min_w
190    for i, t in enumerate(f.get_time()):
191        w = f(t, point_id = k)[0]
192        if w > max_w: max_w = w
193        if w < min_w: min_w = w       
194
195    if max_w-z <= 0.01:  # Find first gauge where max depth > eps (runup)
196        runup_point = g[0]
197
198    if min_w-z <= 0.01:  # Find first gauge where min depth > eps (coastline)
199        coastline = g[0]       
200       
201    max_stage.append(max_w)
202    min_stage.append(min_w)   
203
204
205# Print
206print 'wave height [m]:                    ', amplitude
207runup_height = topography([runup_point], [(north+south)/2.])[0]
208print 'run up height [m]:                  ', runup_height
209
210runup_distance = runup_point-coastline
211print 'run up distance from coastline [m]: ', runup_distance
212
213print 'Coastline (meters form west):       ', coastline
214
215
216
217
218# Stop here
219import sys; sys.exit() 
220
221# Take snapshots and plot
222ion()
223figure(1)
224plot(x_vector, topography(x_vector,(north+south)/2.), 'r-')
225xlabel('x')
226ylabel('Elevation')
227#legend(('Max stage', 'Min stage', 'Elevation'), shadow=True, loc='upper right')
228title('Stage snapshots (t=0, 10, ...) for gauge line')
229grid()
230hold(True)
231
232for i, t in enumerate(f.get_time()):
233    if i % 10 == 0:
234        # Take only some timesteps to avoid clutter
235        stages = []   
236        for k, g in enumerate(gauges):
237            w = f(t, point_id = k)[0]       
238            stages.append(w)
239
240        plot(x_vector, stages, 'b-')
241         
242savefig('snapshots')
243
244
245
246# Store
247filename = 'maxrunup'+str(amplitude)+'.csv'
248fid = open(filename,'w')   
249s = 'Waveheight,Runup distance,Runup height\n'
250fid.write(s)
251
252s = '%.2f,%.2f,%.2f\n' %(amplitude, runup_distance, runup_height)
253fid.write(s)
254
255fid.close()
256
257# Plot max runup etc
258ion()
259figure(1)
260plot(x_vector, max_stage, 'g+',
261     x_vector, min_stage, 'b+',     
262     x_vector, topography(x_vector,(north+south)/2.), 'r-')
263xlabel('x')
264ylabel('stage')
265legend(('Max stage', 'Min stage', 'Elevation'), shadow=True, loc='upper right')
266title('Maximum stage for gauge line')
267grid()
268#axis([33000, 47000, -1000, 3000])
269savefig('max_stage')
270
271close('all')
272   
273
274
Note: See TracBrowser for help on using the repository browser.