source: anuga_validation/solitary_waves/solitary_wave_runup.py @ 4631

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

Refactored limit2007 to tight_slope_limiters (using eclipse)

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 = 0.25 # 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=100.,
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.