source: anuga_core/source/anuga/examples/beach.py @ 3656

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

Work on parallel examples

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