source: anuga_validation/solitary_waves/solitary_wave_runup.py @ 4678

Last change on this file since 4678 was 4678, checked in by ole, 15 years ago

Added validation scripts and data from UQ (Matt Barnes). It does not run yet.
Also committed work on solitary_wave_runup.py

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