Changeset 9238
- Timestamp:
- Jul 3, 2014, 10:42:36 AM (10 years ago)
- Location:
- trunk/anuga_core/validation_tests/case_studies/towradgi
- Files:
-
- 2 added
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/validation_tests/case_studies/towradgi/Compare_results_with_fieldObs.py
r9057 r9238 54 54 pyplot.savefig('Error_peakstage.png') 55 55 56 #pyplot.clf()57 #pyplot.scatter(floodLevels[:,0],floodLevels[:,1])58 #pyplot.gca().set_aspect('equal')59 #for i in range(len(modelled_level)):60 # difString=str(round(floodLevels[i,3]-modelled_level[i],2))61 # pyplot.text(floodLevels[i,0], floodLevels[i,1], difString)62 63 56 64 57 # Make a bunch of GIS outputs … … 80 73 pyplot.figure(figsize=(12,6)) 81 74 pyplot.plot([X.min(),X.max()],[Y.min(),Y.max()],' ') 82 print 'Imshow'83 #import pdb84 #pdb.set_trace()85 75 pyplot.imshow(scipy.flipud(myDepth),extent=[X.min(),X.max(),Y.min(),Y.max()],origin='lower',cmap=pyplot.get_cmap('Greys')) 86 76 pyplot.gca().set_aspect('equal') 87 77 pyplot.colorbar(orientation='horizontal').set_label('Peak Depth in model (m)') 88 78 er1=floodLevels[:,3]-modelled_level 89 #er1=er1*(er1<1.0) + 1.0*(er1>=1.0)90 #er1=er1*(er1> -1.0) - 1.0*(er1<=-1.0)91 79 pyplot.scatter(floodLevels[:,0], floodLevels[:,1], c=er1,s=20,cmap=pyplot.get_cmap('spectral')) 92 pyplot.colorbar().set_label(label='Field observation - Modelled Peak Depth(m)')80 pyplot.colorbar().set_label(label='Field observation - Modelled Peak Stage (m)') 93 81 pyplot.xlim([p.x.min()+p.xllcorner,p.x.max()+p.xllcorner]) 94 82 pyplot.ylim([p.y.min()+p.yllcorner,p.y.max()+p.yllcorner]) -
trunk/anuga_core/validation_tests/case_studies/towradgi/Towradgi_historic_flood.py
r9174 r9238 2 2 Towradgi Creek 17 August 1998 Storm Event Calibration 3 3 Ubuntu Linux 12.04 LTS 64 bit 4 By Petar Milevski, minor edits by Gareth Davies4 By Petar Milevski, some revisions by Gareth Davies 5 5 """ 6 6 #------------------------------------------------------------------------------ … … 43 43 #------------------------------------------------------------------------------ 44 44 45 def read_polygon_dir(weight_dict, directory, filepattern='*.csv'): 46 # get a list of matching filenames in the directory 47 pattern = os.path.join(directory, filepattern) 48 files = glob.glob(pattern) 49 50 # check that the dictionary contains *all* the files 51 52 errors = [] 53 for f in files: 54 try: 55 _ = weight_dict[f] 56 except KeyError: 57 errors.append(f) 58 59 if errors: 60 msg = '' 61 for f in errors: 62 msg = msg + ', ' + f 63 raise KeyError, 'Files not defined in dictionary: %s' % msg[2:] 64 65 # now get the result list 45 def read_polygon_list(poly_list): 46 # Alternative to read_polygon_dir -- allows us to control order of polygons 66 47 result = [] 67 for f in files:68 result.append((read_polygon( f), weight_dict[f]))48 for i in range(len(poly_list)): 49 result.append((read_polygon(poly_list[i][0]) , poly_list[i][1])) 69 50 return result 70 51 … … 73 54 #------------------------------------------------------------------------------ 74 55 75 CatchmentDictionary = { 76 join('Model', 'Bdy', 'CreekBanks.csv'): 16.0, 77 join('Model', 'Bdy', 'FineCatchment.csv'): 36.0, 78 join('Model', 'Bdy', 'Catchment.csv'): 100.0 79 } 80 81 ManningDictionary = { 82 join('Model', 'Mannings', '1.csv'): 0.03, #park 83 join('Model', 'Mannings', '2.csv'): 0.15, 84 join('Model', 'Mannings', '3.csv'): 0.15, 85 join('Model', 'Mannings', '4.csv'): 0.04, 86 join('Model', 'Mannings', '5.csv'): 0.15, 87 join('Model', 'Mannings', '6.csv'): 0.15, 88 join('Model', 'Mannings', '7.csv'): 0.15, 89 join('Model', 'Mannings', '8.csv'): 0.15, 90 join('Model', 'Mannings', '9.csv'): 0.03, #park 91 join('Model', 'Mannings', '10.csv'): 0.15, 92 join('Model', 'Mannings', '11.csv'): 0.15, 93 join('Model', 'Mannings', '12.csv'): 0.15, 94 join('Model', 'Mannings', '13.csv'): 0.025, 95 join('Model', 'Mannings', '14.csv'): 0.15, 96 join('Model', 'Mannings', '15.csv'): 0.15, 97 join('Model', 'Mannings', '16.csv'): 0.15, 98 join('Model', 'Mannings', '17.csv'): 0.15, 99 join('Model', 'Mannings', '18.csv'): 0.045, 100 join('Model', 'Mannings', '18a.csv'): 0.15, 101 join('Model', 'Mannings', '18b.csv'): 0.15, 102 join('Model', 'Mannings', '18c.csv'): 0.15, 103 join('Model', 'Mannings', '18d.csv'): 0.15, 104 join('Model', 'Mannings', '18e.csv'): 0.08, #cokeworks site 105 join('Model', 'Mannings', '19.csv'): 0.15, 106 join('Model', 'Mannings', '20.csv'): 0.15, 107 join('Model', 'Mannings', '21.csv'): 0.15, 108 join('Model', 'Mannings', '22.csv'): 0.15, 109 join('Model', 'Mannings', '23.csv'): 0.15, 110 join('Model', 'Mannings', '24.csv'): 0.05, 111 join('Model', 'Mannings', '25.csv'): 0.15, 112 join('Model', 'Mannings', '26.csv'): 0.15, 113 join('Model', 'Mannings', '27.csv'): 0.15, 114 join('Model', 'Mannings', '28.csv'): 0.15, 115 join('Model', 'Mannings', '29.csv'): 0.15, 116 join('Model', 'Mannings', '30.csv'): 0.15, 117 join('Model', 'Mannings', '31.csv'): 0.15, 118 join('Model', 'Mannings', '32.csv'): 0.15, 119 join('Model', 'Mannings', '33.csv'): 0.15, 120 join('Model', 'Mannings', '34.csv'): 0.15, 121 join('Model', 'Mannings', '35.csv'): 0.15, 122 join('Model', 'Mannings', '36.csv'): 0.05, 123 join('Model', 'Mannings', '37.csv'): 0.15, 124 join('Model', 'Mannings', '38.csv'): 0.15, 125 join('Model', 'Mannings', '39.csv'): 0.15, 126 join('Model', 'Mannings', '40.csv'): 0.15, 127 join('Model', 'Mannings', '41.csv'): 0.15, 128 join('Model', 'Mannings', '42.csv'): 0.15, 129 join('Model', 'Mannings', '43.csv'): 0.15, 130 join('Model', 'Mannings', '44.csv'): 0.15, 131 join('Model', 'Mannings', '45.csv'): 0.15, 132 join('Model', 'Mannings', '46.csv'): 0.15, 133 join('Model', 'Mannings', '47.csv'): 0.15, 134 join('Model', 'Mannings', '48.csv'): 0.15, 135 join('Model', 'Mannings', '49.csv'): 0.15, 136 join('Model', 'Mannings', '50.csv'): 0.15, 137 join('Model', 'Mannings', '51.csv'): 0.15, 138 join('Model', 'Mannings', '52.csv'): 0.15, 139 join('Model', 'Mannings', '53.csv'): 0.15, 140 join('Model', 'Mannings', '54.csv'): 0.15, 141 join('Model', 'Mannings', '55.csv'): 0.15, 142 join('Model', 'Mannings', '56.csv'): 0.15, 143 join('Model', 'Mannings', '57.csv'): 0.15, 144 join('Model', 'Mannings', '58.csv'): 0.15, 145 join('Model', 'Mannings', '59.csv'): 0.08, 146 join('Model', 'Mannings', '60.csv'): 0.15, 147 join('Model', 'Mannings', '61.csv'): 0.08, 148 join('Model', 'Mannings', '62.csv'): 0.15, 149 join('Model', 'Mannings', '63.csv'): 0.08, 150 join('Model', 'Mannings', '64.csv'): 0.15, 151 join('Model', 'Mannings', '65.csv'): 0.15, 152 join('Model', 'Mannings', '66.csv'): 0.15, 153 join('Model', 'Mannings', '67.csv'): 0.15, 154 join('Model', 'Mannings', '68.csv'): 0.15, 155 join('Model', 'Mannings', '69.csv'): 0.15, 156 join('Model', 'Mannings', '70.csv'): 0.15, 157 join('Model', 'Mannings', '71.csv'): 0.05, 158 join('Model', 'Mannings', '72.csv'): 0.15, 159 join('Model', 'Mannings', '73.csv'): 0.15, 160 join('Model', 'Mannings', '74.csv'): 0.15, 161 join('Model', 'Mannings', '75.csv'): 0.15, 162 join('Model', 'Mannings', '76.csv'): 0.15, 163 join('Model', 'Mannings', '77.csv'): 0.07, 164 join('Model', 'Mannings', '78.csv'): 0.15, 165 join('Model', 'Mannings', '79.csv'): 0.15, 166 join('Model', 'Mannings', '80.csv'): 0.15, 167 join('Model', 'Mannings', '81.csv'): 0.15, 168 join('Model', 'Mannings', '82.csv'): 0.15, 169 join('Model', 'Mannings', '83.csv'): 0.15, 170 join('Model', 'Mannings', '84.csv'): 0.15, 171 join('Model', 'Mannings', '85.csv'): 0.15, 172 join('Model', 'Mannings', '86.csv'): 0.15, 173 join('Model', 'Mannings', 'Escarpement.csv'): 0.15, 174 join('Model', 'Mannings', 'Railway.csv'): 0.025, 175 join('Model', 'Creeks', 'creeks.csv'): 0.025 176 } 177 178 basename = join('DEM', 'towradgi') 56 CatchmentList = [ 57 [join('Model', 'Bdy', 'Catchment.csv'), 100.0], 58 [join('Model', 'Bdy', 'FineCatchment.csv'), 36.0], 59 [join('Model', 'Bdy', 'CreekBanks.csv'), 8.0] 60 ] 61 62 ## IMPORTANT -- The ORDER in ManningList matters: When there is overlap, 63 ## priority regions at BOTTOM 64 ## FIXME: This setup can be done with fewer lines of code! 65 channel_manning=0.03 66 ManningList = [ 67 [ join('Model', 'Mannings', '1.csv'),0.04], #park 68 [ join('Model', 'Mannings', '2.csv'),0.15], 69 [ join('Model', 'Mannings', '3.csv'),0.15], 70 [ join('Model', 'Mannings', '4.csv'),0.04], 71 [ join('Model', 'Mannings', '5.csv'),0.15], 72 [ join('Model', 'Mannings', '6.csv'),0.15], 73 [ join('Model', 'Mannings', '7.csv'),0.15], 74 [ join('Model', 'Mannings', '8.csv'),0.15], 75 [ join('Model', 'Mannings', '9.csv'),0.04], #park 76 [ join('Model', 'Mannings', '10.csv'), 0.15], 77 [ join('Model', 'Mannings', '11.csv'), 0.15], 78 [ join('Model', 'Mannings', '12.csv'), 0.15], 79 [ join('Model', 'Mannings', '13.csv'), 0.04], 80 [ join('Model', 'Mannings', '14.csv'), 0.15], 81 [ join('Model', 'Mannings', '15.csv'), 0.15], 82 [ join('Model', 'Mannings', '16.csv'), 0.15], 83 [ join('Model', 'Mannings', '17.csv'), 0.15], 84 [ join('Model', 'Mannings', '18.csv'), 0.045], 85 [ join('Model', 'Mannings', '18a.csv'), 0.15], 86 [ join('Model', 'Mannings', '18b.csv'), 0.15], 87 [ join('Model', 'Mannings', '18c.csv'), 0.15], 88 [ join('Model', 'Mannings', '18d.csv'), 0.15], 89 [ join('Model', 'Mannings', '18e.csv'), 0.08], #cokeworks site 90 [ join('Model', 'Mannings', '19.csv'), 0.15], 91 [ join('Model', 'Mannings', '20.csv'), 0.15], 92 [ join('Model', 'Mannings', '21.csv'), 0.15], 93 [ join('Model', 'Mannings', '22.csv'), 0.15], 94 [ join('Model', 'Mannings', '23.csv'), 0.15], 95 [ join('Model', 'Mannings', '24.csv'), 0.05], 96 [ join('Model', 'Mannings', '25.csv'), 0.15], 97 [ join('Model', 'Mannings', '26.csv'), 0.15], 98 [ join('Model', 'Mannings', '27.csv'), 0.15], 99 [ join('Model', 'Mannings', '28.csv'), 0.15], 100 [ join('Model', 'Mannings', '29.csv'), 0.15], 101 [ join('Model', 'Mannings', '30.csv'), 0.15], 102 [ join('Model', 'Mannings', '31.csv'), 0.15], 103 [ join('Model', 'Mannings', '32.csv'), 0.15], 104 [ join('Model', 'Mannings', '33.csv'), 0.15], 105 [ join('Model', 'Mannings', '34.csv'), 0.15], 106 [ join('Model', 'Mannings', '35.csv'), 0.15], 107 [ join('Model', 'Mannings', '36.csv'), 0.05], 108 [ join('Model', 'Mannings', '37.csv'), 0.15], 109 [ join('Model', 'Mannings', '38.csv'), 0.15], 110 [ join('Model', 'Mannings', '39.csv'), 0.15], 111 [ join('Model', 'Mannings', '40.csv'), 0.15], 112 [ join('Model', 'Mannings', '41.csv'), 0.15], 113 [ join('Model', 'Mannings', '42.csv'), 0.15], 114 [ join('Model', 'Mannings', '43.csv'), 0.15], 115 [ join('Model', 'Mannings', '44.csv'), 0.15], 116 [ join('Model', 'Mannings', '45.csv'), 0.15], 117 [ join('Model', 'Mannings', '46.csv'), 0.15], 118 [ join('Model', 'Mannings', '47.csv'), 0.15], 119 [ join('Model', 'Mannings', '48.csv'), 0.15], 120 [ join('Model', 'Mannings', '49.csv'), 0.15], 121 [ join('Model', 'Mannings', '50.csv'), 0.15], 122 [ join('Model', 'Mannings', '51.csv'), 0.15], 123 [ join('Model', 'Mannings', '52.csv'), 0.15], 124 [ join('Model', 'Mannings', '53.csv'), 0.15], 125 [ join('Model', 'Mannings', '54.csv'), 0.15], 126 [ join('Model', 'Mannings', '55.csv'), 0.15], 127 [ join('Model', 'Mannings', '56.csv'), 0.15], 128 [ join('Model', 'Mannings', '57.csv'), 0.15], 129 [ join('Model', 'Mannings', '58.csv'), 0.15], 130 [ join('Model', 'Mannings', '59.csv'), 0.08], 131 [ join('Model', 'Mannings', '60.csv'), 0.15], 132 [ join('Model', 'Mannings', '61.csv'), 0.08], 133 [ join('Model', 'Mannings', '62.csv'), 0.15], 134 [ join('Model', 'Mannings', '63.csv'), 0.08], 135 [ join('Model', 'Mannings', '64.csv'), 0.15], 136 [ join('Model', 'Mannings', '65.csv'), 0.15], 137 [ join('Model', 'Mannings', '66.csv'), 0.15], 138 [ join('Model', 'Mannings', '67.csv'), 0.15], 139 [ join('Model', 'Mannings', '68.csv'), 0.15], 140 [ join('Model', 'Mannings', '69.csv'), 0.15], 141 [ join('Model', 'Mannings', '70.csv'), 0.15], 142 [ join('Model', 'Mannings', '71.csv'), 0.05], 143 [ join('Model', 'Mannings', '72.csv'), 0.15], 144 [ join('Model', 'Mannings', '73.csv'), 0.15], 145 [ join('Model', 'Mannings', '74.csv'), 0.15], 146 [ join('Model', 'Mannings', '75.csv'), 0.15], 147 [ join('Model', 'Mannings', '76.csv'), 0.15], 148 [ join('Model', 'Mannings', '77.csv'), 0.07], 149 [ join('Model', 'Mannings', '78.csv'), 0.15], 150 [ join('Model', 'Mannings', '79.csv'), 0.15], 151 [ join('Model', 'Mannings', '80.csv'), 0.15], 152 [ join('Model', 'Mannings', '81.csv'), 0.15], 153 [ join('Model', 'Mannings', '82.csv'), 0.15], 154 [ join('Model', 'Mannings', '83.csv'), 0.15], 155 [ join('Model', 'Mannings', '84.csv'), 0.15], 156 [ join('Model', 'Mannings', '85.csv'), 0.15], 157 [ join('Model', 'Mannings', '86.csv'), 0.15], 158 [ join('Model', 'Mannings', 'Escarpement.csv'), 0.15], 159 [ join('Model', 'Mannings', 'Railway.csv'), 0.04], 160 [ join('Model', 'Creeks', 'creeks.csv'), channel_manning] 161 ] 162 163 basename = join('DEM_bridges', 'towradgi') 179 164 outname = join('Towradgi_historic_flood') 180 meshname = join('DEM ','towradgi.tsh')165 meshname = join('DEM_bridges','towradgi.tsh') 181 166 182 167 W=303517 … … 187 172 maximum_triangle_area = 1000 188 173 #minimum_storable_height = 0.05 189 base_friction = 0.0 2174 base_friction = 0.04 190 175 alpha = 0.99 191 176 verbose = True 192 177 193 #model_output_dir='MODEL_OUTPUTS_'+time.strftime('%Y%m%d_%H%M%S', time.localtime())194 178 model_output_dir='MODEL_OUTPUTS' 195 179 try: … … 206 190 207 191 bounding_polygon = [[W, S], [E, S], [E, N], [W, N]] 208 interior_regions = read_polygon_dir(CatchmentDictionary, join('Model', 'Bdy')) 209 192 #interior_regions = read_polygon_dir(CatchmentDictionary, join('Model', 'Bdy')) 193 interior_regions = read_polygon_list(CatchmentList) 194 195 # FIXME: Have these in a shapefile / other file and read them in 196 breaklines=[[[306612.336559443,6193708.75358065], 197 [306604.441364239,6193693.17994946]], 198 [[306977.886673843,6193753.44134088], 199 [306978.027867398,6193710.94208076]], 200 [[306956.001672788,6193750.89985688], 201 [306956.707640564,6193706.14149989]], 202 [[306627.303076293,6193697.45809624], 203 [306620.525785644,6193683.62112783]], 204 [[307236.83565407,6193741.01630802], 205 [307231.682089306,6193721.03741996]], 206 [[307224.975395434,6193742.71063068], 207 [307220.880782334,6193723.36711362]], 208 [[307624.764946969,6193615.98941489], 209 [307617.98765632,6193601.44647871]], 210 [[307613.328268998,6193623.19028621], 211 [307607.751123568,6193610.97704368]]] 212 213 # Make the mesh 210 214 create_mesh_from_regions(bounding_polygon, 211 215 boundary_tags={'south': [0], 'east': [1], 'north': [2], 'west': [3]}, … … 213 217 interior_regions=interior_regions, 214 218 filename=meshname, 219 breaklines=breaklines, 215 220 use_cache=False, 216 221 verbose=True) … … 221 226 222 227 domain = Domain(meshname, use_cache=False, verbose=True) 223 224 from anuga.utilities.argparsing import parse_standard_args 225 alg, cfl = parse_standard_args() 228 # Get arguments 229 args = anuga.get_args() 230 alg = args.alg 231 #cfl = args.cfl 232 verbose = args.verbose 226 233 domain.set_flow_algorithm(alg) 227 234 #domain.set_CFL(cfl) 228 #domain.set_flow_algorithm('DE1')229 #domain.set_store_vertices_smoothly()230 #domain.set_minimum_storable_height(minimum_storable_height) 235 if(not domain.get_using_discontinuous_elevation()): 236 raise Exception, 'This model run relies on a discontinuous elevation solver (because of how topography is set up)' 237 231 238 domain.set_datadir(model_output_dir) 232 239 try: … … 242 249 #------------------------------------------------------------------------------ 243 250 244 friction_list = read_polygon_ dir(ManningDictionary, join('Model', 'Mannings'))251 friction_list = read_polygon_list(ManningList) 245 252 domain.set_quantity('friction', Polygon_function(friction_list, default=base_friction, geo_reference=domain.geo_reference)) 246 253 … … 249 256 250 257 # Decompress the zip file to make a csv for reading 251 #zipfile.ZipFile('DEM/towradgi.zip').extract('towradgi.csv',path='DEM/') 252 # Set the elevation with the new csv file 253 domain.set_quantity('elevation', filename=basename+'.csv', use_cache=True, verbose=True, alpha=alpha) 254 #os.remove('DEM/towradgi.csv') # Clean up csv file 258 zipfile.ZipFile('DEM_bridges/towradgi_cleaner.zip').extract('towradgi.csv',path='DEM_bridges/') 259 260 from anuga.utilities.quantity_setting_functions import make_nearestNeighbour_quantity_function 261 elev_xyz=numpy.genfromtxt(fname=basename+'.csv',delimiter=',') 262 263 # Use nearest-neighbour interpolation of elevation 264 elev_fun_wrapper=make_nearestNeighbour_quantity_function(elev_xyz,domain) 265 domain.set_quantity('elevation', elev_fun_wrapper, location='centroids') 266 267 os.remove('DEM_bridges/towradgi.csv') # Clean up csv file 255 268 else: 256 269 domain = None 257 270 258 271 barrier() 259 if myid == 0 and verbose: print 'DISTRIBUTING DOMAIN' 272 if myid == 0 and verbose: 273 print 'DISTRIBUTING DOMAIN' 274 260 275 domain = distribute(domain) 261 276 barrier() 277 262 278 263 if myid == 0 and verbose: print 'CREATING INLETS' 279 if myid == 0 and verbose: 280 print 'CREATING INLETS' 264 281 """ 282 #FIXME: Include these again 283 265 284 #------------------------------------------------------------------------------ 266 285 # ENTER CULVERT DATA … … 270 289 el0 = numpy.array([[305772.982,6193988.557] , [305772.378,6193987.823]]) 271 290 el1 = numpy.array([[305794.592,6193983.907] , [305793.988,6193983.173]]) 272 273 culvert = Boyd_box_operator(domain, 274 losses=losses, 275 width=0.9, 276 exchange_lines=[el0, el1], 277 height=0.9, 278 apron=3.0, 279 enquiry_gap=10.0, 280 use_momentum_jet=False, 281 use_velocity_head=False, 282 manning=0.013, 283 logging=False, 284 label=join('Runs', '1H'), # this puts the culvert hydraulics output into the same dir as the sww's 285 verbose=False) 291 292 ## Adjust el0, el1 293 #elOffset=0. 294 #el0M=0.5*(el0[0,:]+el0[1,:]) ; el1M=0.5*(el1[0,:]+el1[1,:]); n0=el0M-el1M; n0=n0/((n0*n0).sum())**0.5; 295 #el0 = el0 296 culvert = Boyd_box_operator(domain, 297 losses=losses, 298 width=0.9, 299 exchange_lines=[el0, el1], 300 height=0.9, 301 apron=3.0, 302 enquiry_gap=10.0, 303 use_momentum_jet=True, 304 use_velocity_head=True, 305 manning=0.013, 306 logging=False, 307 label=join('Runs', '1H'), # this puts the culvert hydraulics output into the same dir as the sww's 308 verbose=False) 286 309 287 310 # MeadowStCulvert branch 2 … … 291 314 292 315 culvert = Boyd_box_operator(domain, 293 294 295 296 297 298 299 use_momentum_jet=False,300 use_velocity_head=False,301 302 303 304 316 losses=losses, 317 width=5.4, 318 exchange_lines=[el0, el1], 319 height=0.6, 320 apron=3.0, 321 enquiry_gap=10.0, 322 use_momentum_jet=True, 323 use_velocity_head=True, 324 manning=0.013, 325 logging=False, 326 label=join('Runs', '2H'), # this puts the culvert hydraulics output into the same dir as the sww's 327 verbose=False) 305 328 306 329 # WilliamsStCulvert branch 2 … … 310 333 311 334 culvert = Boyd_box_operator(domain, 312 313 314 315 316 317 318 use_momentum_jet=False,319 use_velocity_head=False,320 321 322 323 335 losses=losses, 336 width=1.2, 337 exchange_lines=[el0, el1], 338 height=1.2, 339 apron=3.0, 340 enquiry_gap=10.0, 341 use_momentum_jet=True, 342 use_velocity_head=True, 343 manning=0.013, 344 logging=False, 345 label=join('Runs', '3H'), # this puts the culvert hydraulics output into the same dir as the sww's 346 verbose=False) 324 347 325 348 # MeadowStCulvert Towradgi branch … … 329 352 330 353 culvert = Boyd_box_operator(domain, 331 332 333 334 335 336 337 use_momentum_jet=False,338 use_velocity_head=False,339 340 341 342 354 losses=losses, 355 width=4.0, 356 exchange_lines=[el0, el1], 357 height=2.2, 358 apron=3.0, 359 enquiry_gap=10.0, 360 use_momentum_jet=True, 361 use_velocity_head=True, 362 manning=0.013, 363 logging=False, 364 label=join('Runs', '4H'), # this puts the culvert hydraulics output into the same dir as the sww's 365 verbose=False) 343 366 344 367 # CollinsStCulverts tarra 5 branch … … 348 371 349 372 culvert = Boyd_box_operator(domain, 350 351 352 353 354 355 356 use_momentum_jet=False,357 use_velocity_head=False,358 359 360 361 373 losses=losses, 374 width=3.6, 375 exchange_lines=[el0, el1], 376 height=0.93, 377 apron=3.0, 378 enquiry_gap=10.0, 379 use_momentum_jet=True, 380 use_velocity_head=True, 381 manning=0.013, 382 logging=False, 383 label=join('Runs', '9H'), # this puts the culvert hydraulics output into the same dir as the sww's 384 verbose=False) 362 385 363 386 # Norther Distributor Culverts branch 5 … … 367 390 368 391 culvert = Boyd_box_operator(domain, 369 370 371 372 373 374 375 use_momentum_jet=False,376 use_velocity_head=False,377 378 379 380 392 losses=losses, 393 width=9.0, 394 exchange_lines=[el0, el1], 395 height=0.85, 396 apron=3.0, 397 enquiry_gap=10.0, 398 use_momentum_jet=True, 399 use_velocity_head=True, 400 manning=0.013, 401 logging=False, 402 label=join('Runs', '10H'), # this puts the culvert hydraulics output into the same dir as the sww's 403 verbose=False) 381 404 382 405 #CokeWorksCulverts branch 5 … … 386 409 387 410 culvert = Boyd_box_operator(domain, 388 389 390 391 392 393 394 use_momentum_jet=False,395 use_velocity_head=False,396 397 398 399 411 losses=losses, 412 width=6.2, 413 exchange_lines=[el0, el1], 414 height=3.0, 415 apron=3.1, 416 enquiry_gap=10.0, 417 use_momentum_jet=True, 418 use_velocity_head=True, 419 manning=0.013, 420 logging=False, 421 label=join('Runs', '11H'), # this puts the culvert hydraulics output into the same dir as the sww's 422 verbose=False) 400 423 401 424 #Northern Distributor Branch 6 Culverts … … 405 428 406 429 culvert = Boyd_box_operator(domain, 407 408 409 410 411 412 413 use_momentum_jet=False,414 use_velocity_head=False,415 416 417 418 430 losses=losses, 431 width=3.6, 432 exchange_lines=[el0, el1], 433 height=1.2, 434 apron=3.1, 435 enquiry_gap=10.0, 436 use_momentum_jet=True, 437 use_velocity_head=True, 438 manning=0.013, 439 logging=False, 440 label=join('Runs', '12H'), # this puts the culvert hydraulics output into the same dir as the sww's 441 verbose=False) 419 442 420 443 #Railway Branch 6 Culverts … … 424 447 425 448 culvert = Boyd_box_operator(domain, 426 427 428 429 430 431 432 use_momentum_jet=False,433 use_velocity_head=False,434 435 436 437 449 losses=losses, 450 width=3.6, 451 exchange_lines=[el0, el1], 452 height=1.2, 453 apron=3.1, 454 enquiry_gap=10.0, 455 use_momentum_jet=True, 456 use_velocity_head=True, 457 manning=0.013, 458 logging=False, 459 label=join('Runs', '13H'), # this puts the culvert hydraulics output into the same dir as the sww's 460 verbose=False) 438 461 439 462 #Colgong Branch6 Culverts … … 443 466 444 467 culvert = Boyd_box_operator(domain, 445 446 447 448 449 450 451 use_momentum_jet=False,452 use_velocity_head=False,453 454 455 456 468 losses=losses, 469 width=2.1, 470 exchange_lines=[el0, el1], 471 height=1.0, 472 apron=3.1, 473 enquiry_gap=10.0, 474 use_momentum_jet=True, 475 use_velocity_head=True, 476 manning=0.013, 477 logging=False, 478 label=join('Runs', '14H'), # this puts the culvert hydraulics output into the same dir as the sww's 479 verbose=False) 457 480 458 481 #Basin Outlet Branch 3 Culverts … … 462 485 463 486 culvert = Boyd_box_operator(domain, 464 465 466 467 468 469 470 use_momentum_jet=False,471 use_velocity_head=False,472 473 474 475 487 losses=losses, 488 width=6.0, 489 exchange_lines=[el0, el1], 490 height=0.86, 491 apron=3.1, 492 enquiry_gap=10.0, 493 use_momentum_jet=True, 494 use_velocity_head=True, 495 manning=0.013, 496 logging=False, 497 label=join('Runs', '15H'), # this puts the culvert hydraulics output into the same dir as the sww's 498 verbose=False) 476 499 477 500 #BellambiRd Branch 3 Culverts … … 481 504 482 505 culvert = Boyd_box_operator(domain, 483 484 485 486 487 488 489 use_momentum_jet=False,490 use_velocity_head=False,491 492 493 494 506 losses=losses, 507 width=1.05, 508 exchange_lines=[el0, el1], 509 height=1.0, 510 apron=3.1, 511 enquiry_gap=10.0, 512 use_momentum_jet=True, 513 use_velocity_head=True, 514 manning=0.013, 515 logging=False, 516 label=join('Runs', '16H'), # this puts the culvert hydraulics output into the same dir as the sww's 517 verbose=False) 495 518 496 519 #MeadowSt Branch 3 Culverts … … 500 523 501 524 culvert = Boyd_box_operator(domain, 502 503 504 505 506 507 508 use_momentum_jet=False,509 use_velocity_head=False,510 511 512 513 525 losses=losses, 526 width=1.5, 527 exchange_lines=[el0, el1], 528 height=1.0, 529 apron=3.1, 530 enquiry_gap=10.0, 531 use_momentum_jet=True, 532 use_velocity_head=True, 533 manning=0.013, 534 logging=False, 535 label=join('Runs', '17H'), # this puts the culvert hydraulics output into the same dir as the sww's 536 verbose=False) 514 537 515 538 #13MeadowSt Branch 3 Culverts … … 519 542 520 543 culvert = Boyd_box_operator(domain, 521 522 523 524 525 526 527 use_momentum_jet=False,528 use_velocity_head=False,529 530 531 532 544 losses=losses, 545 width=1.5, 546 exchange_lines=[el0, el1], 547 height=1.0, 548 apron=3.1, 549 enquiry_gap=10.0, 550 use_momentum_jet=True, 551 use_velocity_head=True, 552 manning=0.013, 553 logging=False, 554 label=join('Runs', '18H'), # this puts the culvert hydraulics output into the same dir as the sww's 555 verbose=False) 533 556 534 557 #41AngelStBranch3Culverts … … 538 561 539 562 culvert = Boyd_box_operator(domain, 540 541 542 543 544 545 546 use_momentum_jet=False,547 use_velocity_head=False,548 549 550 551 563 losses=losses, 564 width=10.0, 565 exchange_lines=[el0, el1], 566 height=0.35, 567 apron=3.1, 568 enquiry_gap=10.0, 569 use_momentum_jet=True, 570 use_velocity_head=True, 571 manning=0.013, 572 logging=False, 573 label=join('Runs', '19H'), # this puts the culvert hydraulics output into the same dir as the sww's 574 verbose=False) 552 575 553 576 #CarrollSt Branch 7 Culverts … … 557 580 558 581 culvert = Boyd_box_operator(domain, 559 560 561 562 563 564 565 use_momentum_jet=False,566 use_velocity_head=False,567 568 569 570 582 losses=losses, 583 width=1.22, 584 exchange_lines=[el0, el1], 585 height=0.3, 586 apron=3.1, 587 enquiry_gap=10.0, 588 use_momentum_jet=True, 589 use_velocity_head=True, 590 manning=0.013, 591 logging=False, 592 label=join('Runs', '20H'), # this puts the culvert hydraulics output into the same dir as the sww's 593 verbose=False) 571 594 572 595 #ParkerRd Branch 7 Culverts … … 576 599 577 600 culvert = Boyd_box_operator(domain, 578 579 580 581 582 583 584 use_momentum_jet=False,585 use_velocity_head=False,586 587 588 589 601 losses=losses, 602 width=3.2, 603 exchange_lines=[el0, el1], 604 height=0.3, 605 apron=3.1, 606 enquiry_gap=10.0, 607 use_momentum_jet=True, 608 use_velocity_head=True, 609 manning=0.013, 610 logging=False, 611 label=join('Runs', '21H'), # this puts the culvert hydraulics output into the same dir as the sww's 612 verbose=False) 590 613 591 614 #LakePde Branch 7 Culverts 592 615 losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0} 593 el0 = numpy.array([[308251.257,6193614.658] , [308248.343,6193614.482]]) 594 el1 = numpy.array([[308233.627,6193598.088] , [308230.713,6193597.912]]) 595 596 culvert = Boyd_box_operator(domain, 597 losses=losses, 598 width=3.0, 599 exchange_lines=[el0, el1], 600 height=0.75, 601 apron=3.1, 602 enquiry_gap=10.0, 603 use_momentum_jet=False, 604 use_velocity_head=False, 605 manning=0.013, 606 logging=False, 607 label=join('Runs', '22H'), # this puts the culvert hydraulics output into the same dir as the sww's 608 verbose=False) 609 616 el0 = numpy.array([[308251.257,6193614.658] , [308248.343,6193618.]]) 617 el1 = numpy.array([[308232.,6193593.] , [308225.,6193596.]]) 618 619 culvert = Boyd_box_operator(domain, 620 losses=losses, 621 width=3.0, 622 exchange_lines=[el0, el1], 623 height=0.75, 624 apron=3.1, 625 enquiry_gap=10.0, 626 use_momentum_jet=True, 627 use_velocity_head=True, 628 manning=0.013, 629 logging=False, 630 label=join('Runs', '22H'), # this puts the culvert hydraulics output into the same dir as the sww's 631 verbose=False) 632 """ 633 634 635 smoothTS=30. # Smoothing timescale for bridges 610 636 # PrincesHwyBridge Towradgi branch 611 losses = {'inlet':0. 5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0}612 el0 = numpy.array([[306608. 128,6193709.170] , [306600.772,6193696.949]])613 el1 = numpy.array([[30662 7.238,6193696.530] , [306619.882,6193684.309]])614 615 culvert = Boyd_box_operator(domain,616 losses=losses,617 width=12.0,618 exchange_lines=[el0, el1],619 height=3.0,620 apron=3.0,621 enquiry_gap=10.0,622 use_momentum_jet=False,623 use_velocity_head=False,624 manning=0.013,625 626 627 637 losses = {'inlet':0.0, 'outlet':0.0, 'bend':0.0, 'grate':0.0, 'pier': 1.0, 'other': 0.0} 638 el0 = numpy.array([[306608.,6193703.] , [306607.,6193700.0]]) 639 el1 = numpy.array([[306624.,6193692.7] , [306623.,6193688.]]) 640 culvert = Boyd_box_operator(domain, 641 losses=losses, 642 width=12.0, 643 exchange_lines=[el0, el1], 644 height=3.0, 645 apron=0.0, 646 enquiry_gap=10.0, 647 smoothing_timescale=smoothTS, 648 use_momentum_jet=True, 649 use_velocity_head=True, 650 manning=channel_manning, 651 logging=False, 652 label=join('Runs', '5H'), # this puts the culvert hydraulics output into the same dir as the sww's 653 verbose=False) 628 654 629 655 # PioneerRdBridge Towradgi branch 630 losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0} 631 el0 = numpy.array([[307613.615,6193631.615] , [307597.940,6193608.357]]) 632 el1 = numpy.array([[307635.865,6193620.848] , [307267.195,6193705.076]]) 633 634 culvert = Boyd_box_operator(domain, 635 losses=losses, 636 width=20.0, 637 exchange_lines=[el0, el1], 638 height=3.0, 639 apron=3.0, 640 enquiry_gap=10.0, 641 use_momentum_jet=False, 642 use_velocity_head=False, 643 manning=0.013, 644 logging=False, 645 label=join('Runs', '8H'), # this puts the culvert hydraulics output into the same dir as the sww's 646 verbose=False) 656 losses = {'inlet':0.0, 'outlet':0.0, 'bend':0.0, 'grate':0.0, 'pier': 1.0, 'other': 0.0} 657 el0 = numpy.array([[307623.,6193610.] , [307622.,6193607.]]) 658 el1 = numpy.array([[307610.,6193619.] , [307609., 6193616.]]) 659 enq0 = numpy.array([[ 307637., 6193588. ]]) 660 enq1 = numpy.array([[ 307596., 6193623. ]]) 661 culvert = Boyd_box_operator(domain, 662 losses=losses, 663 width=20.0, 664 exchange_lines=[el0, el1], 665 #enquiry_points=[enq0, enq1], # This seemed to make stability worse 666 height=3.5, 667 apron=0.0, 668 enquiry_gap=10.0, 669 smoothing_timescale=smoothTS, 670 use_momentum_jet=True, 671 use_velocity_head=True, 672 manning=channel_manning, 673 logging=False, 674 label=join('Runs', '8H'), # this puts the culvert hydraulics output into the same dir as the sww's 675 verbose=False) 647 676 648 677 # NorthernDistributorBridge Towradgi branch 649 losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0} 650 el0 = numpy.array([[306944.945,6193760.601] , [306944.569,6193699.449]]) 651 el1 = numpy.array([[306989.167,6193757.821] , [306987.147,6193699.441]]) 652 653 culvert = Boyd_box_operator(domain, 654 losses=losses, 655 width=60.0, 656 exchange_lines=[el0, el1], 657 height=6.0, 658 apron=3.0, 659 enquiry_gap=10.0, 660 use_momentum_jet=False, 661 use_velocity_head=False, 662 manning=0.013, 663 logging=False, 664 label=join('Runs', '6H'), # this puts the culvert hydraulics output into the same dir as the sww's 665 verbose=False) 666 667 # RailwayBridge Towradgi branch, THIS IS THE BRIDGE THAT BLOCKEDAND DAM BREAK OCCURRED 668 losses = {'inlet':0.5, 'outlet':1.0, 'bend':0.0, 'grate':0.0, 'pier': 0.0, 'other': 0.0} 669 el0 = numpy.array([[307209.474,6193752.810] , [307197.942,6193716.493]]) 670 el1 = numpy.array([[307275.853,6193735.391] , [307228.492,6193726.613]]) 671 672 culvert = Boyd_box_operator(domain, 673 losses=losses, 674 width=3.0, 675 exchange_lines=[el0, el1], 676 height=1.0, 677 apron=3.0, 678 enquiry_gap=10.0, 679 use_momentum_jet=False, 680 use_velocity_head=False, 681 manning=0.013, 682 logging=False, 683 label=join('Runs', '7H'), # this puts the culvert hydraulics output into the same dir as the sww's 684 verbose=False) 685 """ 678 losses = {'inlet':0.0, 'outlet':0.0, 'bend':0.0, 'grate':0.0, 'pier': 1.0, 'other': 0.0} 679 el0 = numpy.array([[306985.,6193749.] , [306985.,6193736.]]) 680 el1 = numpy.array([[306950.,6193745.] , [306950.,6193732.]]) 681 682 #enq0 = numpy.array([[ 306996., 6193750.]]) 683 #enq1 = numpy.array([[ 306931., 6193738.]]) 684 culvert = Boyd_box_operator(domain, 685 losses=losses, 686 width=45.0, 687 exchange_lines=[el0, el1], 688 #enquiry_points=[enq0, enq1], 689 height=6.0, 690 apron=0.0, 691 enquiry_gap=10., #None, #10.0, 692 smoothing_timescale=smoothTS, 693 use_momentum_jet=True, 694 use_velocity_head=True, 695 manning=channel_manning, 696 logging=False, 697 label=join('Runs', '6H'), # this puts the culvert hydraulics output into the same dir as the sww's 698 verbose=False) 699 700 # RailwayBridge Towradgi branch ? (THIS IS THE BRIDGE THAT BLOCKEDAND DAM BREAK OCCURRED) ? 701 losses = {'inlet':0.0, 'outlet':0.0, 'bend':0.0, 'grate':0.0, 'pier': 1.0, 'other': 0.0} 702 el0 = numpy.array([[307236.,6193737.] , [307235.,6193733.]]) 703 el1 = numpy.array([[307223.,6193738.] , [307222.,6193734.]]) 704 culvert = Boyd_box_operator(domain, 705 losses=losses, 706 width=20.0, 707 exchange_lines=[el0, el1], 708 height=8.0, 709 apron=0.0, 710 enquiry_gap=20.0, 711 smoothing_timescale=smoothTS, 712 use_momentum_jet=True, 713 use_velocity_head=True, 714 manning=channel_manning, 715 logging=False, 716 label=join('Runs', '7H'), # this puts the culvert hydraulics output into the same dir as the sww's 717 verbose=False) 718 686 719 687 720 #------------------------------------------------------------------------------ … … 689 722 #------------------------------------------------------------------------------ 690 723 724 # FIXME: Reduce the number of code lines here! 691 725 Catchment_Rain_Polygon100 = read_polygon(join('Forcing', 'Rainfall', 'Gauge', '100.csv')) 692 726 rainfall100 = file_function(join('Forcing', 'Rainfall', 'Hort', '100.tms'), quantities=['rate']) … … 948 982 domain.set_boundary({'west': Bd, 'south': Bd, 'north': Bd, 'east': Bw}) 949 983 950 #----------------------------------------------------------------951 # Change conditions during the run to simulate Pioneer Rd Bridge952 #----------------------------------------------------------------953 #954 #ManningDictionary = {join('Model', 'Creeks', 'creeks.csv'): 1.0}955 #friction_list = read_polygon_dir(ManningDictionary, join('Model', 'Creeks'))956 #957 #ManningDictionary = {join('Model', 'Creeks', 'creeks.csv'): 0.025}958 #friction_list = read_polygon_dir(ManningDictionary, join('Model', 'Creeks'))959 #960 #import pdb961 #pdb.set_trace()962 #963 #for i in range(10):964 # print 'stalling'965 966 984 #------------------------------------------------------------------------------ 967 985 # EVOLVE SYSTEM THROUGH TIME … … 975 993 for t in domain.evolve(yieldstep = 300., finaltime = 83700.): 976 994 #if t == 37800.0: #time when bridge deck starts to get submerged, increase n to act as bridge deck, handrail and blockage effects 977 # print 'Reset friction...' 978 # domain.set_quantity('friction',Polygon_function(friction_list, default=base_friction, geo_reference=domain.geo_reference)) 995 ## Try to block all culverts / bridges, as described in the flood study 979 996 #if t == 44100.0: #time when water level drops below bridge deck, bring n back down to existing conditions 980 997 # print 'Reset friction...' 981 # domain.set_quantity('friction', Polygon_function(friction_list, default=base_friction, geo_reference=domain.geo_reference))982 998 if myid == 0: 983 999 domain.write_time()
Note: See TracChangeset
for help on using the changeset viewer.