# source:anuga_core/source/anuga/examples/beach.py@3563

Last change on this file since 3563 was 3563, checked in by ole, 17 years ago

Moved shallow water out from the old pyvolution directory.
All tests pass, most examples and okushiri works again.

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