1 | """Simple water flow example using ANUGA |
---|
2 | |
---|
3 | Water driven up a linear slope and time varying boundary, |
---|
4 | similar to a beach environment |
---|
5 | """ |
---|
6 | |
---|
7 | |
---|
8 | #------------------------------------------------------------------------------ |
---|
9 | # Import necessary modules |
---|
10 | #------------------------------------------------------------------------------ |
---|
11 | |
---|
12 | from pmesh.mesh_interface import create_mesh_from_regions |
---|
13 | from pyvolution.shallow_water import Domain |
---|
14 | from pyvolution.shallow_water import Reflective_boundary |
---|
15 | from pyvolution.shallow_water import Dirichlet_boundary |
---|
16 | from pyvolution.shallow_water import Time_boundary |
---|
17 | from pyvolution.shallow_water import Transmissive_boundary |
---|
18 | |
---|
19 | |
---|
20 | #------------------------------------------------------------------------------ |
---|
21 | # Setup computational domain |
---|
22 | #------------------------------------------------------------------------------ |
---|
23 | |
---|
24 | waveheight = 1800.0 |
---|
25 | depth_east_edge = -4000. |
---|
26 | timeend = 1000.0 |
---|
27 | west = 0. |
---|
28 | east = 80000. |
---|
29 | south = 0. |
---|
30 | north = 5000. |
---|
31 | polygon = [[east,north],[west,north],[west,south],[east,south]] |
---|
32 | meshname = 'test.msh' |
---|
33 | |
---|
34 | create_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 | |
---|
43 | domain = Domain(meshname,use_cache=True,verbose = True) |
---|
44 | domain.set_name('test') # Output to test.sww |
---|
45 | domain.set_datadir('.') # Use current directory for output |
---|
46 | domain.set_quantities_to_be_stored(['stage',# Store all conserved quantities |
---|
47 | 'xmomentum', |
---|
48 | 'ymomentum']) |
---|
49 | |
---|
50 | |
---|
51 | #------------------------------------------------------------------------------ |
---|
52 | # Setup initial conditions |
---|
53 | #------------------------------------------------------------------------------ |
---|
54 | |
---|
55 | def 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 | |
---|
60 | domain.set_quantity('elevation', topography) # Use function for elevation |
---|
61 | domain.set_quantity('friction', 0. ) # Constant friction |
---|
62 | domain.set_quantity('stage', 0.0) # Constant initial stage |
---|
63 | |
---|
64 | #------------------------------------------------------------------------------ |
---|
65 | # Setup boundary conditions |
---|
66 | #------------------------------------------------------------------------------ |
---|
67 | |
---|
68 | from math import sin, pi |
---|
69 | Br = Reflective_boundary(domain) # Solid reflective wall |
---|
70 | Bt = Transmissive_boundary(domain) # Continue all values on boundary |
---|
71 | Bd = Dirichlet_boundary([0.,0.,0.]) # Constant boundary values |
---|
72 | Bw = 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 |
---|
76 | domain.set_boundary({'left': Br, 'right': Bw, 'top': Br, 'bottom': Br}) |
---|
77 | |
---|
78 | |
---|
79 | #------------------------------------------------------------------------------ |
---|
80 | # Evolve system through time |
---|
81 | #------------------------------------------------------------------------------ |
---|
82 | |
---|
83 | stagestep = [] |
---|
84 | for t in domain.evolve(yieldstep = 10, finaltime = timeend): |
---|
85 | domain.write_time() |
---|
86 | |
---|
87 | #----------------------------------------------------------------------------- |
---|
88 | # Interrogate solution |
---|
89 | #----------------------------------------------------------------------------- |
---|
90 | |
---|
91 | # Gauge line |
---|
92 | def 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 | |
---|
103 | gauges, gaugex, gaugey = gauge_line(west,east,north,south) |
---|
104 | |
---|
105 | from pyvolution.util import file_function |
---|
106 | |
---|
107 | f = file_function('test.sww', |
---|
108 | quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'], |
---|
109 | interpolation_points = gauges, |
---|
110 | verbose = True, |
---|
111 | use_cache = True) |
---|
112 | |
---|
113 | from pylab import plot, xlabel, ylabel, title, ion, close, savefig, figure, axis |
---|
114 | ion() |
---|
115 | |
---|
116 | maxw = [] |
---|
117 | minw = [] |
---|
118 | maxd = [] |
---|
119 | count = 0 |
---|
120 | for 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 | |
---|
142 | filename = 'maxrunup'+str(waveheight)+str(depth_east_edge)+'.csv' |
---|
143 | fid = open(filename,'w') |
---|
144 | s = 'Depth at eastern boundary,Waveheight,Runup distance,Runup height\n' |
---|
145 | fid.write(s) |
---|
146 | runup_distance = posx |
---|
147 | runup_height = topography(runup_distance,(north+south)/2.) |
---|
148 | print 'run up height: ', runup_height |
---|
149 | if runup_distance < (east+west)/2.: |
---|
150 | runup_distance_from_coast = (east+west)/2.-runup_distance |
---|
151 | else: |
---|
152 | runup_distance_from_coast = runup_distance-(east+west)/2. |
---|
153 | print 'run up distance from coastline: ', runup_distance_from_coast |
---|
154 | s = '%.2f,%.2f,%.2f,%.2f\n' %(depth_east_edge,waveheight,runup_distance_from_coast,runup_height) |
---|
155 | fid.write(s) |
---|
156 | |
---|
157 | figure(1) |
---|
158 | plot(gaugex,maxw,'g+',gaugex,topography(gaugex,(north+south)/2.),'r-') |
---|
159 | xlabel('x') |
---|
160 | ylabel('stage') |
---|
161 | title('Maximum stage for gauge line') |
---|
162 | #axis([33000, 47000, -1000, 3000]) |
---|
163 | savefig('max_stage') |
---|
164 | |
---|
165 | close('all') |
---|