Changeset 1351
- Timestamp:
- May 9, 2005, 12:23:58 PM (20 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/ga/storm_surge/Merimbula/run_merimbula_lake.py
r1295 r1351 4 4 Christopher Zoppou, Stephen Roberts, Ole Nielsen, Duncan Gray 5 5 Geoscience Australia, ANU 6 6 7 7 Specific methods pertaining to the 2D shallow water equation 8 8 are imported from shallow_water … … 24 24 Dirichlet_boundary, Wind_stress 25 25 from pmesh2domain import pmesh_to_domain_instance 26 from util import file_function, Polygon_function, read_polygon 27 from Numeric import zeros, Float 26 from util import file_function, Polygon_function, read_polygon, inside_polygon 27 from Numeric import zeros, Float, asarray 28 from least_squares import Interpolation 29 import time 28 30 29 31 #------- 30 32 # Domain 31 filename = 'merimbula_interpolated .tsh'32 filename = 'merimbula_10 785_1.tsh'33 filename = 'merimbula_interpolated_bathymetry.tsh' 34 filename = 'merimbula_10834_bridge_refined_bathymetry.tsh' 33 35 print 'Creating domain from', filename 34 36 domain = pmesh_to_domain_instance(filename, Domain) … … 36 38 37 39 #------------------------------------------ 38 # Reduction operation for get_vertex_values 40 # Reduction operation for get_vertex_values 39 41 from util import mean 40 domain.reduction = mean 42 domain.reduction = mean 41 43 42 #-----------------43 # Scale weed zones44 def weed_zone(points):45 #print points46 n = len(points)47 print n48 z = []49 for i in range(n):50 #print i51 z.append([0,0])52 z[i][0] = 755471.4 + (points[i][0] + 3250.)/1.12553 z[i][1] = 5910260.0 + (points[i][1] + 1337.)/1.1254 #print z55 return z56 44 57 #------------- 58 # Set friction 59 # Sand bed 60 w = 0.035 61 # Saltmarsh 62 g = 0.060 63 # Paddle weed 64 y = 0.025 65 # Eel grass 66 r = 0.035 67 # Mangroves 68 c = 0.065 69 # Strap weed 70 b = 0.040 71 #----------------------------------------------- 72 # Set the whole region to a constant value 73 weed_zoneall = weed_zone([[-2269,-1337],[1894,-1339],[1894,2946],[-2669,2946]]) 74 75 #--------------------------------------- 76 # Read friction polygon boundaries 77 weed_zone47 = weed_zone(read_polygon('weed_zone.047')) 78 weed_zone2 = weed_zone(read_polygon('weed_zone.002')) 79 weed_zone12 = weed_zone(read_polygon('weed_zone.012')) 80 weed_zone35 = weed_zone(read_polygon('weed_zone.035')) 81 weed_zone8 = weed_zone(read_polygon('weed_zone.008')) 82 weed_zone10 = weed_zone(read_polygon('weed_zone.010')) 83 weed_zone13 = weed_zone(read_polygon('weed_zone.013')) 84 weed_zone15 = weed_zone(read_polygon('weed_zone.015')) 85 weed_zone19 = weed_zone(read_polygon('weed_zone.019')) 86 weed_zone18 = weed_zone(read_polygon('weed_zone.018')) 87 weed_zone24 = weed_zone(read_polygon('weed_zone.024')) 88 weed_zone26 = weed_zone(read_polygon('weed_zone.026')) 89 weed_zone27 = weed_zone(read_polygon('weed_zone.027')) 90 weed_zone32 = weed_zone(read_polygon('weed_zone.032')) 91 weed_zone31 = weed_zone(read_polygon('weed_zone.031')) 92 weed_zone33 = weed_zone(read_polygon('weed_zone.033')) 93 weed_zone34 = weed_zone(read_polygon('weed_zone.034')) 94 weed_zone36 = weed_zone(read_polygon('weed_zone.036')) 95 weed_zone37 = weed_zone(read_polygon('weed_zone.037')) 96 weed_zone38 = weed_zone(read_polygon('weed_zone.038')) 97 weed_zone40 = weed_zone(read_polygon('weed_zone.040')) 98 weed_zone41 = weed_zone(read_polygon('weed_zone.041')) 99 weed_zone42 = weed_zone(read_polygon('weed_zone.042')) 100 weed_zone43 = weed_zone(read_polygon('weed_zone.043')) 101 weed_zone44 = weed_zone(read_polygon('weed_zone.044')) 102 weed_zone45 = weed_zone(read_polygon('weed_zone.045')) 103 weed_zone46 = weed_zone(read_polygon('weed_zone.046')) 104 weed_zone1 = weed_zone(read_polygon('weed_zone.001')) 105 weed_zone20 = weed_zone(read_polygon('weed_zone.020')) 106 weed_zone21 = weed_zone(read_polygon('weed_zone.021')) 107 weed_zone22 = weed_zone(read_polygon('weed_zone.022')) 108 weed_zone23 = weed_zone(read_polygon('weed_zone.023')) 109 weed_zone25 = weed_zone(read_polygon('weed_zone.025')) 110 weed_zone16 = weed_zone(read_polygon('weed_zone.016')) 111 weed_zone17 = weed_zone(read_polygon('weed_zone.017')) 112 weed_zone3 = weed_zone(read_polygon('weed_zone.003')) 113 weed_zone6 = weed_zone(read_polygon('weed_zone.006')) 114 weed_zone7 = weed_zone(read_polygon('weed_zone.007')) 115 weed_zone9 = weed_zone(read_polygon('weed_zone.009')) 116 weed_zone4 = weed_zone(read_polygon('weed_zone.004')) 117 weed_zone39 = weed_zone(read_polygon('weed_zone.039')) 118 weed_zone28 = weed_zone(read_polygon('weed_zone.028')) 119 weed_zone29 = weed_zone(read_polygon('weed_zone.029')) 120 weed_zone30 = weed_zone(read_polygon('weed_zone.030')) 121 weed_zone5 = weed_zone(read_polygon('weed_zone.005')) 122 weed_zone11 = weed_zone(read_polygon('weed_zone.011')) 123 weed_zone14 = weed_zone(read_polygon('weed_zone.014')) 124 125 print weed_zoneall 126 127 domain.set_quantity('friction',Polygon_function( 128 [ (weed_zoneall, w), (weed_zone47, b), (weed_zone2, y), (weed_zone12, y) ] )) 129 130 ## [(weed_zone35, 'y')],[(weed_zone8, 'r')],[(weed_zone10, 'g')],[(weed_zone13, 'g')], \ 131 ## [(weed_zone15, 'g')],[(weed_zone19, 'c')],[(weed_zone18, 'g')],[(weed_zone24, 'c')], \ 132 ## [(weed_zone26, 'r')],[(weed_zone27, 'g')],[(weed_zone32, 'c')],[(weed_zone31, 'g')], \ 133 ## [(weed_zone33, 'r')],[(weed_zone34, 'g')],[(weed_zone36, 'r')],[(weed_zone37, 'g')], \ 134 ## [(weed_zone38, 'g')],[(weed_zone40, 'r')],[(weed_zone41, 'b')],[(weed_zone42, 'r')], \ 135 ## [(weed_zone43, 'r')],[(weed_zone44, 'r')],[(weed_zone45, 'b')],[(weed_zone46, 'r')], \ 136 ## [(weed_zone1, 'w')],[(weed_zone20, 'w')],[(weed_zone21, 'w')],[(weed_zone22, 'w')], \ 137 ## [(weed_zone23, 'w')],[(weed_zone25, 'w')],[(weed_zone16, 'w')],[(weed_zone17, 'w')], \ 138 ## [(weed_zone3, 'w')],[(weed_zone6, 'w')],[(weed_zone7, 'w')],[(weed_zone9, 'w')], \ 139 ## [(weed_zone4, 'b')],[(weed_zone39, 'b')],[(weed_zone28, 'r')],[(weed_zone29, 'b')], \ 140 ## [(weed_zone30, 'b')],[(weed_zone5, 'b')],[(weed_zone11, 'b')],[(weed_zone14, 'b')]]) 141 142 ## [[(weed_zoneall,'w')],[(weed_zone47, 'b')],[(weed_zone2, 'y')],[(weed_zone12, 'y')], \ 143 ## [(weed_zone35, 'y')],[(weed_zone8, 'r')],[(weed_zone10, 'g')],[(weed_zone13, 'g')], \ 144 ## [(weed_zone15, 'g')],[(weed_zone19, 'c')],[(weed_zone18, 'g')],[(weed_zone24, 'c')], \ 145 ## [(weed_zone26, 'r')],[(weed_zone27, 'g')],[(weed_zone32, 'c')],[(weed_zone31, 'g')], \ 146 ## [(weed_zone33, 'r')],[(weed_zone34, 'g')],[(weed_zone36, 'r')],[(weed_zone37, 'g')], \ 147 ## [(weed_zone38, 'g')],[(weed_zone40, 'r')],[(weed_zone41, 'b')],[(weed_zone42, 'r')], \ 148 ## [(weed_zone43, 'r')],[(weed_zone44, 'r')],[(weed_zone45, 'b')],[(weed_zone46, 'r')], \ 149 ## [(weed_zone1, 'w')],[(weed_zone20, 'w')],[(weed_zone21, 'w')],[(weed_zone22, 'w')], \ 150 ## [(weed_zone23, 'w')],[(weed_zone25, 'w')],[(weed_zone16, 'w')],[(weed_zone17, 'w')], \ 151 ## [(weed_zone3, 'w')],[(weed_zone6, 'w')],[(weed_zone7, 'w')],[(weed_zone9, 'w')], \ 152 ## [(weed_zone4, 'b')],[(weed_zone39, 'b')],[(weed_zone28, 'r')],[(weed_zone29, 'b')], \ 153 ## [(weed_zone30, 'b')],[(weed_zone5, 'b')],[(weed_zone11, 'b')],[(weed_zone14, 'b')]])) 154 45 domain.set_quantity('friction',0.03) 155 46 #-------------------- 156 47 # Boundary conditions 48 157 49 #--------------------------------------- 158 50 # Tidal cycle recorded at Eden as open … … 160 52 print 'Open sea boundary condition from ',filename 161 53 Bf = File_boundary(filename, domain) 54 162 55 #-------------------------------------- 163 56 # All other boundaries are reflective … … 175 68 #-------------------------------- 176 69 # Initial water surface elevation 177 domain.set_quantity('stage', 0.0)178 70 domain.set_quantity('stage', -50.0) 71 179 72 #---------------------------------------------------------- 180 73 # Decide which quantities are to be stored at each timestep … … 185 78 domain.store = True #Store for visualisation purposes 186 79 domain.format = 'sww' #Native netcdf visualisation format 187 domain.filename = 'Merimbula_2003 '80 domain.filename = 'Merimbula_2003_4days_dry' 188 81 189 82 #---------------------- 190 83 # Set order of accuracy 191 84 domain.default_order = 1 192 domain.visualise = True193 domain.visualise_color_stage = True194 domain.visualise_stage_range = 1.0195 85 domain.smooth = True 86 print domain.set_inscribed_circle() 196 87 197 88 #--------- 198 #Evolution 199 yieldstep = 10 200 finaltime = 1000 89 # Evolution 90 t0 = time.time() 91 yieldstep = 900 92 finaltime = 588000*3 201 93 for t in domain.evolve(yieldstep = yieldstep, finaltime = finaltime): 202 94 domain.write_time() 203 204 print 'Done' 205 206 95 96 print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset
for help on using the changeset viewer.
