source: documentation/user_manual/examples/runuptest.py @ 3563

Last change on this file since 3563 was 3563, checked in by ole, 16 years ago

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.

File size: 5.5 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# Import necessary modules
10#------------------------------------------------------------------------------
11
12from anuga.pmesh.mesh_interface import create_mesh_from_regions
13from anuga.shallow_water import Domain
14from anuga.shallow_water import Reflective_boundary
15from anuga.shallow_water import Dirichlet_boundary
16from anuga.shallow_water import Time_boundary
17from anuga.shallow_water import Transmissive_boundary
18
19
20#------------------------------------------------------------------------------
21# Setup computational domain
22#------------------------------------------------------------------------------
23
24waveheight = 1800.0
25depth_east_edge = -4000.
26timeend = 1000.0
27west = 0.
28east = 80000.
29south = 0.
30north = 5000.
31polygon = [[east,north],[west,north],[west,south],[east,south]]
32meshname = 'test.msh'
33
34create_mesh_from_regions(polygon,
35                         boundary_tags={'top': [0],
36                                        'left': [1],
37                                        'bottom': [2],
38                                        'right': [3]},
39                         maximum_triangle_area=200000,
40                         filename=meshname)
41                         #interior_regions=interior_regions)
42
43domain = Domain(meshname,use_cache=True,verbose = True)
44domain.set_name('test')                     # Output to test.sww
45domain.set_datadir('.')                     # Use current directory for output
46domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities
47                                    'xmomentum',
48                                    'ymomentum'])   
49
50
51#------------------------------------------------------------------------------
52# Setup initial conditions
53#------------------------------------------------------------------------------
54
55def topography(x,y):
56    slope = depth_east_edge/((east-west)/2.)
57    c = -(west+east)/2.*slope
58    return slope*x+c                         # linear bed slope
59
60domain.set_quantity('elevation', topography) # Use function for elevation
61domain.set_quantity('friction', 0. )         # Constant friction
62domain.set_quantity('stage', 0.0)            # Constant initial stage
63
64#------------------------------------------------------------------------------
65# Setup boundary conditions
66#------------------------------------------------------------------------------
67
68from math import sin, pi
69Br = Reflective_boundary(domain)      # Solid reflective wall
70Bt = Transmissive_boundary(domain)    # Continue all values on boundary
71Bd = Dirichlet_boundary([0.,0.,0.])   # Constant boundary values
72Bw = Time_boundary(domain=domain,     # Time dependent boundary 
73                   f=lambda t: [((10<t<20)*waveheight), 0.0, 0.0])
74
75# Associate boundary tags with boundary objects
76domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br})
77
78
79#------------------------------------------------------------------------------
80# Evolve system through time
81#------------------------------------------------------------------------------
82
83stagestep = []
84for t in domain.evolve(yieldstep = 10, finaltime = timeend):
85    domain.write_time()
86
87#-----------------------------------------------------------------------------
88# Interrogate solution
89#-----------------------------------------------------------------------------
90
91# Gauge line
92def gauge_line(west,east,north,south):
93    from Numeric import arange
94    gaugex = arange(west,(west+east)/2.,1000.)
95    gauges = []
96    gaugey = []
97    for x in gaugex:
98        y = (north+south)/2.
99        gaugey.append(y)
100        gauges.append([x,y])
101    return gauges, gaugex, gaugey
102
103gauges, gaugex, gaugey = gauge_line(west,east,north,south)
104
105from anuga.abstract_2d_finite_volumes.util import file_function
106
107f = file_function('test.sww',
108                  quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'],
109                  interpolation_points = gauges,
110                  verbose = True,
111                  use_cache = True)
112
113from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis
114ion()
115
116maxw = []
117minw = []
118maxd = []
119count = 0
120for k, g in enumerate(gauges):
121    stage = []
122    elevation = []
123    depths = []
124    bed = topography(g[0],g[1])
125    for i, t in enumerate(f.get_time()):
126        w = f(t, point_id = k)[0]
127        elev = f(t, point_id = k)[1]
128        depth = w-elev
129        stage.append(w)
130        elevation.append(elev)
131        depths.append(elev)
132
133    if max(stage)-bed <= 0.2:
134        count+=1
135        posx=g[0]
136        loc = k 
137    maxw.append(max(stage))
138    minw.append(min(stage))
139    maxd.append(max(depths))
140   
141
142filename = 'maxrunup'+str(waveheight)+str(depth_east_edge)+'.csv'
143fid = open(filename,'w')   
144s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n'
145fid.write(s)
146runup_distance = posx
147runup_height = topography(runup_distance,(north+south)/2.)
148print 'run up height:   ', runup_height
149if runup_distance < (east+west)/2.:
150    runup_distance_from_coast = (east+west)/2.-runup_distance
151else:
152    runup_distance_from_coast = runup_distance-(east+west)/2.
153print 'run up distance from coastline: ', runup_distance_from_coast
154s = '%.2f,%.2f,%.2f,%.2f\n' %(depth_east_edge,waveheight,runup_distance_from_coast,runup_height)
155fid.write(s)
156
157figure(1)
158plot(gaugex,maxw,'g+',gaugex,topography(gaugex,(north+south)/2.),'r-')
159xlabel('x')
160ylabel('stage')
161title('Maximum stage for gauge line')
162#axis([33000, 47000, -1000, 3000])
163savefig('max_stage')
164
165close('all')
Note: See TracBrowser for help on using the repository browser.