source: inundation/ga/storm_surge/examples/beach.py @ 1452

Last change on this file since 1452 was 901, checked in by ole, 20 years ago

First experiment with new limiter (near bed)

File size: 5.4 KB
Line 
1"""Example of the inundationmodel.
2
3A wave of water is washed up ontp a hypothetical beach.
4
5To run:
6
7python beach.py
8"""
9
10######################
11# Module imports
12#
13
14import sys
15from os import sep, path
16sys.path.append('..'+sep+'pyvolution')
17
18from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
19     Transmissive_boundary, Time_boundary, Wind_stress
20
21from pmesh2domain import pmesh_to_domain_instance
22from util import read_polygon, Polygon_function
23from math import pi
24from Numeric import choose, greater, ones, sin, exp
25import time
26
27######################
28# Domain
29name = 'beach'
30print 'Creating domain from %s.tsh' %name
31domain = pmesh_to_domain_instance(name + '.tsh', Domain)
32
33#HACK to remove internal boundary (until we get a solution)
34#Is this still necessary? NO
35#for b in domain.boundary.keys():
36#    if domain.boundary[b] == '':
37#        del domain.boundary[b]
38#domain = Domain(domain.coordinates, domain.triangles, domain.boundary)
39
40domain.store = True
41domain.set_name(name + '6')
42domain.default_order = 2
43print "Output being written to " + domain.get_datadir() + sep + \
44      domain.filename + "_size%d." %len(domain) + domain.format
45
46
47def bathymetry(x, y):
48    cut = 75
49   
50    res = choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) ))
51    res == (100-y)/50 + 1
52    #res += (100-y)/130  #Lift up southern bathymetry
53    #res -= y/1500  #Pull down northern bathymetry
54   
55    return res
56   
57
58def topography(x, y):
59
60    import RandomArray
61    #z = 4*sin(pi*x/50) + (100-y)/50 + 1 #Ridge
62   
63    z = -4.0*x/25 + 8  + (100-y)/50 #Beach
64   
65    #z += sin(pi*x/5)/4
66    #z += RandomArray.normal(0, 1, x.shape)/4
67   
68    z += 6*exp( -((x-30)/10)**2 ) * exp( -((y-50)/8)**2 )     #Mound
69    z += 4*exp( -((x-30)/4)**2 ) * exp( -((y-26)/10)**2 )     #extra ridge   
70    z += 4*exp( -((x-30)/5)**2 ) * exp( -((y-10)/8)**2 )     #extra ridge
71    z -= 4*exp( -((x-15)/6)**2 ) * exp( -((y-20)/12)**2 )     #depression
72
73    z += 1.2*exp( -((x-88)/7)**2 ) + exp( -((y-20)/25)**2 )     #Seafloor
74    return z
75
76def riverbed(x, y):
77    return (y-100)/70 - 4.0*x/25 + 8
78   
79
80shoreline = [[40,0], [100,0], [100,100], [65,100], 
81             [55,90], [55,70], [56,50], [50,25], [40,0]]       
82
83land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0],
84        [0,0], [0,100]]
85
86
87water = [[55,0], [100,0], [100,100], [55,100]]
88all = [[0,0], [0,100], [100,100], [100,0]]
89
90                       
91             
92building1 = [[45,80], [49,78], [52,83], [46,83]]             
93building2 = [[35,75], [40,75], [40,80], [35,80]]             
94building3 = [[42,63], [46,61], [48,65], [44,67]]             
95building4 = [[28,56], [28,51], [33,51], [33,56]]             
96building5 = [[10,70], [10,65], [15,65], [15,70]]             
97building6 = [[10,50], [10,45], [15,45], [15,50]]             
98
99river = [[20,100], [18,90], [20,80], [20,60], [15,40], [11,20], [2,0], [10,0],
100         [14,10], [20,30], [24,45], [27,80], [27,85], [35,90], [39,100]]
101             
102print 'Set elevation'   
103t0 = time.time()
104domain.set_quantity('elevation',
105    Polygon_function( [(all, topography),
106                      #(shoreline, bathymetry), (land, topography),
107                      (building1, 7), (building2, 8), 
108                      (building3, 7), (building4, 13),
109                      (building5, 10), (building6, 11)]))
110                      #(river, riverbed)]))
111                     
112print 'That took %.2f seconds' %(time.time()-t0)
113                     
114print 'Set stage'                     
115domain.set_quantity('stage', 
116                    Polygon_function( [(water, -1.5), 
117                    (land, -10)] )) 
118                   
119print 'Set friction'                                       
120domain.set_quantity('friction', 0.08)
121
122print domain.get_extent()
123
124
125
126#Add lateral wind gusts bearing 135 degrees
127def gust(t,x,y): 
128    from math import sin, pi
129    from Numeric import zeros, ones, Float
130
131    N = len(x)
132
133    tt = sin(2*pi*t/50)
134
135    if tt > 0.98:
136        return 26000*tt*ones(N, Float)
137    else:
138        return zeros(N, Float)
139   
140domain.forcing_terms.append(Wind_stress(gust, 135))
141
142
143#Add lateral wind gusts bearing 90 degrees
144def gust2(t,x,y): 
145    from math import sin, pi
146    from Numeric import zeros, ones, Float
147
148    N = len(x)
149
150    tt = sin(2*pi*t/100)
151
152    if tt > 0.95:
153        return 22000*tt*ones(N, Float)
154    else:
155        return zeros(N, Float)
156   
157domain.forcing_terms.append(Wind_stress(gust2, 90))
158
159#Add lateral wind gusts bearing 255 degrees
160def gust3(t,x,y): 
161    from math import sin, pi
162    from Numeric import zeros, ones, Float
163
164    N = len(x)
165
166    tt = sin(2*pi*(t-30)/55)
167
168    if tt > 0.96:
169        return 24000*tt*ones(N, Float)
170    else:
171        return zeros(N, Float)
172   
173domain.forcing_terms.append(Wind_stress(gust3, 255))
174
175
176######################
177# Boundary conditions
178
179print 'Boundaries'
180Br = Reflective_boundary(domain)
181Bo = Transmissive_boundary(domain)
182
183#Constant outflow
184Bd = Dirichlet_boundary([-10, 0.0, 0.0])
185Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0])
186
187print 'Available boundary tags are', domain.get_boundary_tags()
188
189#Set boundary conditions
190tags = {}
191tags['ocean'] = Bt
192tags['wall'] = Br
193tags['wall1'] = Br
194tags['outflow'] = Bd
195tags['exterior'] = Br
196tags['external'] = Br
197#tags['land'] = Bo 
198tags['land'] = Bd 
199tags['westbank'] = None    #Riverbank
200tags['eastbank'] = None    #Riverbank
201tags['eastbankN'] = None    #Riverbank
202domain.set_boundary(tags)
203
204
205domain.check_integrity()
206
207######################
208#Evolution
209t0 = time.time()
210for t in domain.evolve(yieldstep = 0.2, finaltime = 300):
211    domain.write_time()
212print 'Simulation took %.2f seconds' %(time.time()-t0)
213   
Note: See TracBrowser for help on using the repository browser.