source: trunk/anuga_core/documentation/user_manual/examples/runuptest.py @ 7808

Last change on this file since 7808 was 7808, checked in by hudson, 12 years ago

Examples work with new API.

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