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

Last change on this file since 3740 was 3700, checked in by ole, 17 years ago

Improved runup script for Rosh

File size: 7.1 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
6"""
7
8
9#------------------------------------------------------------------------------
10# Import necessary modules
11#------------------------------------------------------------------------------
12
13import sys
14
15from anuga.pmesh.mesh_interface import create_mesh_from_regions
16from abstract_2d_finite_volumes.mesh_factory import rectangular_cross
17from anuga.config import g
18from anuga.shallow_water import Domain
19from anuga.shallow_water import Reflective_boundary
20from anuga.shallow_water import Dirichlet_boundary
21from anuga.shallow_water import Time_boundary
22from anuga.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,\
25     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,
33                    # and h=0 (coast) at 300m
34highest_point = 6   # Highest elevation (m)
35sea_level = 0       # Mean sea level
36min_elevation = -20 # Lowest elevation (elevation of offshore flat part)
37offshore_depth = sea_level-min_elevation # offshore water depth
38
39amplitude = 0.5     # Solitary wave height H
40normalized_amplitude = amplitude/offshore_depth
41simulation_name = 'runup_convergence'   
42coastline_x = -highest_point/slope
43
44# Basin dimensions (m)
45west = 0          # left boundary
46east = 1500       # right boundary
47south = 0         # lower boundary
48north = 100       # upper boundary
49
50
51#------------------------------------------------------------------------------
52# Setup computational domain all units in meters
53#------------------------------------------------------------------------------
54
55# Structured mesh
56dx = 30           # Resolution: Length of subdivisions on x axis (length)
57dy = 30           # Resolution: Length of subdivisions on y axis (width)
58
59length = east-west
60width = north-south
61points, vertices, boundary = rectangular_cross(length/dx, width/dy,
62                                               len1=length, len2=width,
63                                               origin = (west, south)) 
64
65domain = Domain(points, vertices, boundary) # Create domain
66
67
68# Unstructured mesh
69polygon = [[east,north],[west,north],[west,south],[east,south]]
70interior_polygon = [[400,north-10],[west+10,north-10],
71                    [west+10,south+10],[400,south+10]]
72meshname = simulation_name + '.msh'
73create_mesh_from_regions(polygon,
74                         boundary_tags={'top': [0], 'left': [1],
75                                        'bottom': [2], 'right': [3]},
76                         maximum_triangle_area=dx*dy/4,
77                         filename=meshname,
78                         interior_regions=[[interior_polygon,dx*dy/32]])
79
80domain = Domain(meshname, use_cache=True, verbose = True)
81domain.set_minimum_storable_height(0.01)
82
83domain.set_name(simulation_name)
84
85
86#------------------------------------------------------------------------------
87# Setup initial conditions
88#------------------------------------------------------------------------------
89
90#def topography(x,y):
91#    return slope*x+highest_point  # Return linear bed slope (vector)
92
93def topography(x,y):
94    """Two part topography - slope and flat part
95    """
96
97    from Numeric import zeros, Float
98
99    z = zeros(len(x), Float) # Allocate space for return vector
100    for i in range(len(x)):
101
102        z[i] = slope*x[i]+highest_point # Linear bed slope bathymetry
103
104        if z[i] < min_elevation:        # Limit depth
105            z[i] = min_elevation
106
107    return z       
108       
109
110       
111
112domain.set_quantity('elevation', topography) # Use function for elevation
113domain.set_quantity('friction', 0.0 )        # Constant friction
114domain.set_quantity('stage', sea_level)      # Constant initial stage
115
116
117#------------------------------------------------------------------------------
118# Setup boundary conditions
119#------------------------------------------------------------------------------
120
121from math import sin, pi, cosh, sqrt
122Br = Reflective_boundary(domain)      # Solid reflective wall
123Bd = Dirichlet_boundary([0.,0.,0.])   # Constant boundary values
124
125
126def waveform(t): 
127    return sea_level +\
128           amplitude/cosh(((t-50)/offshore_depth)*(0.75*g*amplitude)**0.5)**2
129
130# Time dependent boundary for stage, where momentum is set automatically
131Bts = Transmissive_Momentum_Set_Stage_boundary(domain, waveform)
132
133# Associate boundary tags with boundary objects
134domain.set_boundary({'left': Br, 'right': Bts, 'top': Br, 'bottom': Br})
135
136
137# Find initial runup location and height (coastline)
138w0 = domain.get_maximum_inundation_elevation()
139x0, y0 = domain.get_maximum_inundation_location()
140print
141print 'Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w0, x0, y0)
142
143# Sanity check
144w_i = domain.get_quantity('stage').get_values(interpolation_points=[[x0,y0]])
145print 'Interpolated elevation at (x,y)=(%.2f, %.2f) is %.2f' %(x0, y0, w_i) 
146
147
148#------------------------------------------------------------------------------
149# Evolve system through time
150#------------------------------------------------------------------------------
151
152w_max = w0
153for t in domain.evolve(yieldstep = 1, finaltime = 300):
154    domain.write_time()
155
156    w = domain.get_maximum_inundation_elevation()
157    x, y = domain.get_maximum_inundation_location()
158    print '  Coastline elevation = %.2f at (x,y)=(%.2f, %.2f)' %(w, x, y)
159    print   
160
161    if w > w_max:
162        w_max = w
163        x_max = x
164        y_max = y
165
166
167print '**********************************************'
168print 'Max coastline elevation = %.2f at (%.2f, %.2f)' %(w_max, x_max, y_max)
169print 'Run up distance = %.2f' %sqrt( (x_max-x0)**2 + (y_max-y0)**2 )
170print '**********************************************'
171
172
173#-----------------------------------------------------------------------------
174# Interrogate further
175#---------------------------------------------------------------
176
177# Generate time series of one "gauge" situated at right hand boundary
178from anuga.abstract_2d_finite_volumes.util import sww2timeseries
179production_dirs = {'.': 'test'}
180swwfiles = {}
181for label_id in production_dirs.keys():
182   
183    swwfile = simulation_name + '.sww'
184    swwfiles[swwfile] = label_id
185   
186texname, elev_output = sww2timeseries(swwfiles,
187                                      'boundary_gauge.xya',
188                                      production_dirs,
189                                      report = False,
190                                      reportname = 'test',
191                                      plot_quantity = ['stage', 'speed'],
192                                      surface = False,
193                                      time_min = None,
194                                      time_max = None,
195                                      title_on = True,
196                                      verbose = True)
197
198
199
200
Note: See TracBrowser for help on using the repository browser.