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