Changeset 672 for inundation/ga/storm_surge/examples/beach.py
- Timestamp:
- Dec 6, 2004, 1:02:40 AM (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/examples/beach.py
r668 r672 23 23 from util import read_polygon, Polygon_function 24 24 from math import pi 25 from Numeric import choose, greater, ones, sin 26 25 from Numeric import choose, greater, ones, sin, exp 26 import time 27 27 28 28 ###################### … … 33 33 34 34 #HACK to remove internal boundary (until we get a solution) 35 #Is this still necessary? 35 36 for b in domain.boundary.keys(): 36 37 if domain.boundary[b] == '': … … 46 47 47 48 def bathymetry(x, y): 48 #cut = 5549 #return -(x - cut)/550 51 49 cut = 75 52 return choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) )) 53 ####return choose( greater(x, cut), (-20*ones(x.shape), ones(x.shape))) 50 51 res = choose( greater(x, cut), ( -(x - 55)/5, -4*ones(x.shape) )) 52 res += (100-y)/50 - 0.75 53 return res 54 54 55 55 56 56 def topography(x, y): 57 return 4*sin(pi*x/50) + 1 57 58 import RandomArray 59 #z = 4*sin(pi*x/50) + (100-y)/50 + 1 #Ridge 60 61 z = -4.0*x/25 + 8 + (100-y)/50 #Beach 62 63 #z += sin(pi*x/5)/4 64 #z += RandomArray.normal(0, 1, x.shape)/4 65 66 z += 6*exp( -((x-30)/10)**2 ) * exp( -((y-50)/8)**2 ) #Mound 67 z += 4*exp( -((x-30)/4)**2 ) * exp( -((y-26)/10)**2 ) #extra ridge 68 z += 4*exp( -((x-30)/5)**2 ) * exp( -((y-10)/8)**2 ) #extra ridge 69 70 return z 58 71 59 72 def riverbed(x, y): 60 #return (y-100)/10 + 5*x/100 61 return (y-100)/50 + 4.4 73 return (y-100)/70 - 4.0*x/25 + 8 62 74 63 75 … … 66 78 67 79 land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0], 68 [0,0], [0,100]] 80 [0,0], [0,100]] 81 69 82 70 building1 = [[45,80], [ 50,78], [52,83], [45,83]]83 building1 = [[45,80], [49,78], [52,83], [46,83]] 71 84 building2 = [[35,75], [40,75], [40,80], [35,80]] 72 building3 = [[42,63], [47,63], [50,68], [37,68]] 73 building4 = [[28,70], [28,60], [33,60], [33,70]] 74 building5 = [[10,70], [10,60], [15,60], [15,70]] 75 building6 = [[10,60], [10,50], [15,50], [15,60]] 76 77 78 85 building3 = [[42,63], [46,61], [48,65], [44,67]] 86 building4 = [[28,56], [28,51], [33,51], [33,56]] 87 building5 = [[10,70], [10,65], [15,65], [15,70]] 88 building6 = [[10,50], [10,45], [15,45], [15,50]] 79 89 80 90 river = [[20,100], [18,90], [20,80], [20,60], [15,40], [11,20], [2,0], [10,0], … … 82 92 83 93 print 'Set elevation' 94 t0 = time.time() 84 95 domain.set_quantity('elevation', 85 96 Polygon_function( [(shoreline, bathymetry), (land, topography), 86 (building1, 3), (building2, 9),87 (building3, 8), (building4, 7),88 (building5, 8), (building6, 6),97 (building1, 6), (building2, 8), 98 (building3, 6), (building4, 13), 99 (building5, 10), (building6, 11), 89 100 (river, riverbed)])) 101 102 print 'That took %.2f seconds' %(time.time()-t0) 90 103 91 104 print 'Set level' … … 95 108 96 109 print 'Set friction' 97 domain.set_quantity('friction', 0.0 3)110 domain.set_quantity('friction', 0.075) 98 111 99 112 print domain.get_extent() … … 111 124 112 125 if tt > 0.7: 113 return 10000*tt*ones(N, Float)126 return 30000*tt*ones(N, Float) 114 127 else: 115 128 return zeros(N, Float) … … 125 138 Bo = Transmissive_boundary(domain) 126 139 127 #Constant inflow128 Bd = Dirichlet_boundary([ 0.0, -1.0, 0.0])129 Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0])140 #Constant outflow 141 Bd = Dirichlet_boundary([-10, 0.0, 0.0]) 142 Bt = Time_boundary(domain, lambda t: [ 2.8*(1+sin(2*pi*t/100)), 0.0, 0.0]) 130 143 131 144 print 'Available boundary tags are', domain.get_boundary_tags() … … 135 148 tags['ocean'] = Bt 136 149 tags['wall'] = Br 137 tags['outflow'] = Bo 150 tags['wall1'] = Br 151 tags['outflow'] = Bd 138 152 tags['exterior'] = Br 139 tags['external'] = Br 153 tags['external'] = Br 154 tags['land'] = Bo 155 tags['westbank'] = None #Riverbank 156 tags['eastbank'] = None #Riverbank 157 tags['eastbankN'] = None #Riverbank 140 158 domain.set_boundary(tags) 141 159 … … 145 163 ###################### 146 164 #Evolution 147 for t in domain.evolve(yieldstep = 0.2, finaltime = 500):165 for t in domain.evolve(yieldstep = 0.2, finaltime = 300): 148 166 domain.write_time() 149 167
Note: See TracChangeset
for help on using the changeset viewer.