source: inundation/examples/beach.py @ 3381

Last change on this file since 3381 was 2932, checked in by ole, 19 years ago

Modified beach example to better match the new allowance of non-zero speeds in water shallower than 1mm added in changeset:2561

File size: 5.2 KB
RevLine 
[659]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
[1977]15#sys.path.append('..'+sep+'pyvolution')
[659]16
[2866]17from pyvolution.shallow_water import Domain, Reflective_boundary,\
18     Dirichlet_boundary,\
[668]19     Transmissive_boundary, Time_boundary, Wind_stress
[659]20
[2533]21from pyvolution.pmesh2domain import pmesh_to_domain_instance
22from utilities.polygon import read_polygon, Polygon_function
[667]23from math import pi
[672]24from Numeric import choose, greater, ones, sin, exp
25import time
[659]26
27######################
28# Domain
29name = 'beach'
30print 'Creating domain from %s.tsh' %name
[2866]31domain = Domain(mesh_filename = name + '.tsh',
32                use_cache=True, verbose=True)
[659]33
34domain.store = True
[2568]35#domain.minimum_allowed_height = 0.0
[694]36domain.set_name(name + '6')
[673]37domain.default_order = 2
[901]38print "Output being written to " + domain.get_datadir() + sep + \
[2932]39      domain.filename + domain.format
[659]40
41
[2568]42def bathymetry(x,y):
[667]43    cut = 75
[659]44   
[688]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   
[672]50    return res
51   
[659]52
[667]53def topography(x, y):
[659]54
[672]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   
[688]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
[672]69    return z
70
[667]71def riverbed(x, y):
[672]72    return (y-100)/70 - 4.0*x/25 + 8
[667]73   
[659]74
[667]75shoreline = [[40,0], [100,0], [100,100], [65,100], 
76             [55,90], [55,70], [56,50], [50,25], [40,0]]       
[659]77
[667]78land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0],
[688]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
[672]85                       
[667]86             
[672]87building1 = [[45,80], [49,78], [52,83], [46,83]]             
[667]88building2 = [[35,75], [40,75], [40,80], [35,80]]             
[672]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]]             
[667]93
94river = [[20,100], [18,90], [20,80], [20,60], [15,40], [11,20], [2,0], [10,0],
[668]95         [14,10], [20,30], [24,45], [27,80], [27,85], [35,90], [39,100]]
[667]96             
97print 'Set elevation'   
[672]98t0 = time.time()
[667]99domain.set_quantity('elevation',
[688]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)]))
[667]106                     
[672]107print 'That took %.2f seconds' %(time.time()-t0)
108                     
[774]109print 'Set stage'                     
110domain.set_quantity('stage', 
[688]111                    Polygon_function( [(water, -1.5), 
[2932]112                                       (land, -10)] )) 
[667]113                   
114print 'Set friction'                                       
[2932]115domain.set_quantity('friction', 0.03)
[659]116
117print domain.get_extent()
118
[2582]119print domain.statistics()
[668]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
[688]131    if tt > 0.98:
[2932]132        return 100*tt*ones(N, Float)
[668]133    else:
134        return zeros(N, Float)
135   
[2932]136#domain.forcing_terms.append(Wind_stress(gust, 135))
[668]137
138
[688]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
[694]146    tt = sin(2*pi*t/100)
[688]147
148    if tt > 0.95:
[2932]149        return 100*tt*ones(N, Float)
[688]150    else:
151        return zeros(N, Float)
152   
153domain.forcing_terms.append(Wind_stress(gust2, 90))
154
[694]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
[688]159
[694]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   
[2932]169#domain.forcing_terms.append(Wind_stress(gust3, 255))
[694]170
171
[659]172######################
173# Boundary conditions
174
175print 'Boundaries'
176Br = Reflective_boundary(domain)
[668]177Bo = Transmissive_boundary(domain)
[659]178
[672]179#Constant outflow
180Bd = Dirichlet_boundary([-10, 0.0, 0.0])
[690]181Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0])
[659]182
183print 'Available boundary tags are', domain.get_boundary_tags()
184
185#Set boundary conditions
186tags = {}
[667]187tags['ocean'] = Bt
[659]188tags['wall'] = Br
[672]189tags['wall1'] = Br
190tags['outflow'] = Bd
[659]191tags['exterior'] = Br
[672]192tags['external'] = Br
[673]193#tags['land'] = Bo 
194tags['land'] = Bd 
[672]195tags['westbank'] = None    #Riverbank
196tags['eastbank'] = None    #Riverbank
197tags['eastbankN'] = None    #Riverbank
[659]198domain.set_boundary(tags)
199
200
201domain.check_integrity()
202
203######################
204#Evolution
[688]205t0 = time.time()
[672]206for t in domain.evolve(yieldstep = 0.2, finaltime = 300):
[659]207    domain.write_time()
[688]208print 'Simulation took %.2f seconds' %(time.time()-t0)
[659]209   
Note: See TracBrowser for help on using the repository browser.