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

Last change on this file since 3383 was 3258, checked in by sexton, 19 years ago

fixup

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 pmesh.mesh_interface import create_mesh_from_regions
13from pyvolution.shallow_water import Domain
14from pyvolution.shallow_water import Reflective_boundary
15from pyvolution.shallow_water import Dirichlet_boundary
16from pyvolution.shallow_water import Time_boundary
17from pyvolution.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 pyvolution.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.