[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 | |
---|
[2866] | 17 | from pyvolution.shallow_water import Domain, Reflective_boundary,\ |
---|
| 18 | Dirichlet_boundary,\ |
---|
[668] | 19 | Transmissive_boundary, Time_boundary, Wind_stress |
---|
[659] | 20 | |
---|
[2533] | 21 | from pyvolution.pmesh2domain import pmesh_to_domain_instance |
---|
| 22 | from utilities.polygon import read_polygon, Polygon_function |
---|
[667] | 23 | from math import pi |
---|
[672] | 24 | from Numeric import choose, greater, ones, sin, exp |
---|
| 25 | import time |
---|
[659] | 26 | |
---|
| 27 | ###################### |
---|
| 28 | # Domain |
---|
| 29 | name = 'beach' |
---|
| 30 | print 'Creating domain from %s.tsh' %name |
---|
[2866] | 31 | domain = Domain(mesh_filename = name + '.tsh', |
---|
| 32 | use_cache=True, verbose=True) |
---|
[659] | 33 | |
---|
| 34 | domain.store = True |
---|
[2568] | 35 | #domain.minimum_allowed_height = 0.0 |
---|
[694] | 36 | domain.set_name(name + '6') |
---|
[673] | 37 | domain.default_order = 2 |
---|
[901] | 38 | print "Output being written to " + domain.get_datadir() + sep + \ |
---|
[2932] | 39 | domain.filename + domain.format |
---|
[659] | 40 | |
---|
| 41 | |
---|
[2568] | 42 | def 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] | 53 | def 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] | 71 | def riverbed(x, y): |
---|
[672] | 72 | return (y-100)/70 - 4.0*x/25 + 8 |
---|
[667] | 73 | |
---|
[659] | 74 | |
---|
[667] | 75 | shoreline = [[40,0], [100,0], [100,100], [65,100], |
---|
| 76 | [55,90], [55,70], [56,50], [50,25], [40,0]] |
---|
[659] | 77 | |
---|
[667] | 78 | land = [[65,100], [55,90], [55,70], [56,50], [50,25], [40,0], |
---|
[688] | 79 | [0,0], [0,100]] |
---|
| 80 | |
---|
| 81 | |
---|
| 82 | water = [[55,0], [100,0], [100,100], [55,100]] |
---|
| 83 | all = [[0,0], [0,100], [100,100], [100,0]] |
---|
| 84 | |
---|
[672] | 85 | |
---|
[667] | 86 | |
---|
[672] | 87 | building1 = [[45,80], [49,78], [52,83], [46,83]] |
---|
[667] | 88 | building2 = [[35,75], [40,75], [40,80], [35,80]] |
---|
[672] | 89 | building3 = [[42,63], [46,61], [48,65], [44,67]] |
---|
| 90 | building4 = [[28,56], [28,51], [33,51], [33,56]] |
---|
| 91 | building5 = [[10,70], [10,65], [15,65], [15,70]] |
---|
| 92 | building6 = [[10,50], [10,45], [15,45], [15,50]] |
---|
[667] | 93 | |
---|
| 94 | river = [[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 | |
---|
| 97 | print 'Set elevation' |
---|
[672] | 98 | t0 = time.time() |
---|
[667] | 99 | domain.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] | 107 | print 'That took %.2f seconds' %(time.time()-t0) |
---|
| 108 | |
---|
[774] | 109 | print 'Set stage' |
---|
| 110 | domain.set_quantity('stage', |
---|
[688] | 111 | Polygon_function( [(water, -1.5), |
---|
[2932] | 112 | (land, -10)] )) |
---|
[667] | 113 | |
---|
| 114 | print 'Set friction' |
---|
[2932] | 115 | domain.set_quantity('friction', 0.03) |
---|
[659] | 116 | |
---|
| 117 | print domain.get_extent() |
---|
| 118 | |
---|
[2582] | 119 | print domain.statistics() |
---|
[668] | 120 | |
---|
| 121 | |
---|
| 122 | #Add lateral wind gusts bearing 135 degrees |
---|
| 123 | def 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 |
---|
| 140 | def 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 | |
---|
| 153 | domain.forcing_terms.append(Wind_stress(gust2, 90)) |
---|
| 154 | |
---|
[694] | 155 | #Add lateral wind gusts bearing 255 degrees |
---|
| 156 | def 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 | |
---|
| 175 | print 'Boundaries' |
---|
| 176 | Br = Reflective_boundary(domain) |
---|
[668] | 177 | Bo = Transmissive_boundary(domain) |
---|
[659] | 178 | |
---|
[672] | 179 | #Constant outflow |
---|
| 180 | Bd = Dirichlet_boundary([-10, 0.0, 0.0]) |
---|
[690] | 181 | Bt = Time_boundary(domain, lambda t: [ 3.0*(1+sin(2*pi*t/100)), 0.0, 0.0]) |
---|
[659] | 182 | |
---|
| 183 | print 'Available boundary tags are', domain.get_boundary_tags() |
---|
| 184 | |
---|
| 185 | #Set boundary conditions |
---|
| 186 | tags = {} |
---|
[667] | 187 | tags['ocean'] = Bt |
---|
[659] | 188 | tags['wall'] = Br |
---|
[672] | 189 | tags['wall1'] = Br |
---|
| 190 | tags['outflow'] = Bd |
---|
[659] | 191 | tags['exterior'] = Br |
---|
[672] | 192 | tags['external'] = Br |
---|
[673] | 193 | #tags['land'] = Bo |
---|
| 194 | tags['land'] = Bd |
---|
[672] | 195 | tags['westbank'] = None #Riverbank |
---|
| 196 | tags['eastbank'] = None #Riverbank |
---|
| 197 | tags['eastbankN'] = None #Riverbank |
---|
[659] | 198 | domain.set_boundary(tags) |
---|
| 199 | |
---|
| 200 | |
---|
| 201 | domain.check_integrity() |
---|
| 202 | |
---|
| 203 | ###################### |
---|
| 204 | #Evolution |
---|
[688] | 205 | t0 = time.time() |
---|
[672] | 206 | for t in domain.evolve(yieldstep = 0.2, finaltime = 300): |
---|
[659] | 207 | domain.write_time() |
---|
[688] | 208 | print 'Simulation took %.2f seconds' %(time.time()-t0) |
---|
[659] | 209 | |
---|