source: inundation/examples/beach.py @ 2568

Last change on this file since 2568 was 2568, checked in by ole, 18 years ago

Played with beach example

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