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

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

Renamed pyvolution to abstract_2d_finite_volumes. This is
one step towards pulling pyvolution apart. More to follow.
All unit tests pass and most examples fixed up.




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