Changeset 5503
- Timestamp:
- Jul 15, 2008, 5:00:20 PM (16 years ago)
- Location:
- anuga_work/development/Hinwood_2008
- Files:
-
- 3 edited
-
run_dam.py (modified) (3 diffs)
-
slope.py (modified) (15 diffs)
-
test_slope.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/Hinwood_2008/run_dam.py
r5494 r5503 61 61 friction=0.012, # planed wood. http://www.lmnoeng.com/manningn.htm 62 62 outputdir_name=None, 63 run_type=0): 63 run_type=0, 64 end_tag = '_limiterD'): 64 65 65 66 66 67 basename = 'zz_' + metadata_dic['scenario_id'] 67 end_tag = '_limiterD'68 68 69 if run_type == 1: 69 70 outputdir_name += '_test' + end_tag … … 261 262 # 4 is 0.01 262 263 # 5 is 0.001 263 run_type = 4264 run_type = 5 264 265 #for run_data in [scenarios[5]]: 265 266 #scenarios = scenarios[2:] … … 269 270 run_data, 270 271 run_type = run_type, 271 outputdir_name=run_data['scenario_id']) 272 outputdir_name=run_data['scenario_id'], 273 end_tag='_limiterD') 272 274 #gauges_for_slope(pro_instance.outputdir,[run_data]) -
anuga_work/development/Hinwood_2008/slope.py
r5495 r5503 10 10 11 11 from Numeric import arange, array, zeros, Float, where, greater, less, \ 12 compress, argmin, choose 12 compress, argmin, choose, searchsorted 13 13 14 14 from anuga.fit_interpolate.interpolate import interpolate_sww2csv … … 69 69 70 70 for i in arange(n_slope_locations): 71 slope_locations[i] = (locations[i+1] + locations[i+1])/2.72 71 delta_locations[i] = (locations[i+1] - locations[i]) 72 slope_locations[i] = locations[i] + 0.5*delta_locations[i] 73 73 74 74 for j in arange(n_time): … … 118 118 119 119 # Add the contour line levels to the colorbar 120 cbar.ax.set_ylabel(' verbosity coefficient')120 cbar.ax.set_ylabel('stage slope') 121 121 #cbar.add_lines(CS2) 122 122 … … 146 146 147 147 # Can't seem to reshape this info once it is in the function 148 CS = contourf(x_data, times, z_data, [-1,0.6,0.8,1,2,4], 149 colors = ('black', 'r', 'g', 'b','r'), 150 #cmap=cm.bone, 148 #CS = contourf(x_data, times, z_data, [-1,0.6,0.8,1,2,4], 149 # colors = ('black', 'r', 'g', 'b','r'), 150 # #cmap=cm.bone, 151 # origin=origin) 152 CS = contourf(x_data, times, z_data, 10, 153 #colors = ('black', 'r', 'g', 'b','r'), 154 cmap=cm.bone, 151 155 origin=origin) 152 156 … … 168 172 169 173 # Add the contour line levels to the colorbar 170 cbar.ax.set_ylabel(' verbosity coefficient')174 cbar.ax.set_ylabel('Froude Number') 171 175 #cbar.add_lines(CS2) 172 176 … … 190 194 anuga_break_times.append( \ 191 195 break_time - run_data['ANUGA_start_time']) 192 plot_title = "Stage slope " + id193 196 stage_file = pro_instance.outputdir + "fslope_stage_" + end 197 plot_title = "Stage slope " + id + "\n file:" + \ 198 outputdir_name + "_slope_stage" + plot_type 199 print "Creating ", stage_file 194 200 save_as = pro_instance.plots_dir + sep + \ 195 201 outputdir_name + "_slope_stage" + plot_type 196 202 times, locations, slopes = load_slopes(stage_file) 197 times, slopes = get_band(anuga_break_times[0]-TIME_BORDER,198 100, times, slopes, 0)199 locations, slopes = get_band(200 run_data['break_xs'][0]- 2*LOCATION_BORDER,201 100, locations, slopes, -1)203 #times, slopes = get_band(anuga_break_times[0]-TIME_BORDER, 204 # 100, times, slopes, 0) 205 #locations, slopes = get_band( 206 # min(run_data['break_xs'])- 2*LOCATION_BORDER, 207 # 100, locations, slopes, -1) 202 208 graph_contours(times, locations, slopes, 203 209 plot_title=plot_title, … … 221 227 anuga_break_times.append( \ 222 228 break_time - run_data['ANUGA_start_time']) 223 plot_title = "Froude Number" + id 224 229 plot_title = "Froude Number" + id + "\n file:" + \ 230 outputdir_name + "_froude" + plot_type 225 231 froude_file = pro_instance.outputdir + "fslope_froude_" + end 232 print "Creating ", froude_file 226 233 save_as = pro_instance.plots_dir + sep + \ 227 234 outputdir_name + "_froude" + plot_type … … 229 236 dtimes, sensors = get_band(anuga_break_times[0]-TIME_BORDER, 230 237 100, dtimes, sensors, 0) 231 locations, sensors = get_band(run_data['break_xs'][0]-LOCATION_BORDER, 232 100, locations, sensors, -1) 238 locations, sensors = get_band( 239 min(run_data['break_xs'])-LOCATION_BORDER, 240 100, locations, sensors, -1) 233 241 #print "dtimes", dtimes 234 242 #print "sensors", sensors … … 241 249 is_interactive=is_interactive) 242 250 243 244 def auto_find_min_slopes(outputdir_tag, scenarios): 251 def find_froude(times_froude, locations_froude, froudes_array, 252 times, locations): 253 if len(times) == 0: 254 return [] 255 time_indexes = searchsorted(times_froude, times) 256 location_indexes = searchsorted(locations_froude, locations) 257 258 259 assert len(time_indexes) == len(location_indexes) 260 261 froudes = [] 262 for time_i, loc_i, time, location in map(None, time_indexes, 263 location_indexes, 264 times, locations): 265 # the time values should be the same 266 assert times_froude[time_i] == time 267 268 # The distance value should be half way between the froude locations 269 midpoint = locations_froude[loc_i-1] + \ 270 (locations_froude[loc_i]-locations_froude[loc_i-1])*0.5 271 #print "location", location 272 #print "midpoint", midpoint 273 assert location == midpoint 274 froude = froudes_array[time_i, loc_i-1] + \ 275 (froudes_array[time_i, loc_i]- \ 276 froudes_array[time_i, loc_i-1])*0.5 277 froudes.append(froude) 278 279 return froudes 280 281 def auto_find_min_slopes(slope_tag, outputdir_tag, scenarios): 245 282 """ 246 283 … … 258 295 break_time - run_data['ANUGA_start_time']) 259 296 260 stage_file = pro_instance.outputdir + "fslope_stage_" + end 297 stage_file = pro_instance.outputdir + slope_tag + "slope_stage_" + end 298 froude_file = pro_instance.outputdir + slope_tag + "slope_froude_" + \ 299 end 261 300 262 301 times, slope_locations, slopes = load_slopes(stage_file) 302 #print "slope_locations", slope_locations 303 times_froude, locations_froude, froudes_a = load_sensors(froude_file) 304 #print "locations_froude", locations_froude 263 305 waves = find_min_slopes(times, slope_locations, slopes, 264 run_data['break_times'],265 306 anuga_break_times, 266 307 run_data['band_offset']) 308 267 309 # write the wave info here 268 # return the keys actually. 269 for wave_key in waves.iterkeys(): 270 wave_file = stage_file[:-3] + '_'+ wave_key + ".csv" 310 # and find the froude values 311 for i, wave in enumerate(waves): 312 313 id = "wave_" + str(i) 314 wave_file = stage_file[:-4] + '_'+ id + ".csv" 271 315 print "wave_file", wave_file 316 froudes = find_froude(times_froude, locations_froude, 317 froudes_a, wave[TIME_STR], 318 slope_locations) 272 319 wave_writer = writer(file(wave_file, "wb")) 273 wave_writer.writerow(["x location", "min slope", "Time" ])320 wave_writer.writerow(["x location", "min slope", "Time", "Froude"]) 274 321 wave_writer.writerows(map(None, 275 322 slope_locations, 276 waves[wave_key][SLOPE_STR], 277 waves[wave_key][TIME_STR])) 323 wave[SLOPE_STR], 324 wave[TIME_STR], 325 froudes)) 326 278 327 279 #wave_writer.close() 328 def auto_plot_froude_slopes(slope_tag, outputdir_tag, scenarios): 329 """ 330 331 """ 332 333 plot_type = ".pdf" 334 335 for run_data in scenarios: 336 id = run_data['scenario_id'] 337 outputdir_name = id + outputdir_tag 338 pro_instance = project.Project(['data','flumes','Hinwood_2008'], 339 outputdir_name=outputdir_name) 340 341 assert len(run_data['break_times']) == len(run_data['break_xs']) 342 assert len(run_data['break_times']) == len(run_data['break_type']) 343 344 end = id + ".csv" 345 346 anuga_break_times = [] 347 for break_time in run_data['break_times']: 348 anuga_break_times.append( \ 349 break_time - run_data['ANUGA_start_time']) 280 350 351 #run_data['break_type'] = (run_data['break_type'][0]) 352 for i in range(len(run_data['break_type'])): 281 353 282 283 284 def gauges_for_slope(outputdir_tag, scenarios): 354 wave_tag = "wave_" + str(i) 355 stage_file = pro_instance.outputdir + slope_tag + \ 356 "slope_stage_" + end 357 wave_file = stage_file[:-4] + '_'+ wave_tag + ".csv" 358 save_as = pro_instance.plots_dir + sep + \ 359 outputdir_name + "_" + wave_tag + plot_type 360 print "wave_file", wave_file 361 break_type = run_data['break_type'][i] 362 plot_title = id + ' Wave: ' + str(i) + \ 363 " Break Type: " + break_type + '\n' + \ 364 "File: " + wave_file[34:] # not good! 365 plot_foude_slope_stage(wave_file, 366 anuga_break_times[i], 367 run_data['break_xs'][i], 368 plot_title=plot_title, 369 break_type=break_type, 370 save_as=save_as, 371 is_interactive=False) 372 373 374 375 def gauges_for_slope(slope_tag, outputdir_tag, scenarios): 285 376 """ 286 377 This is used to create a stage file, using gauges relivent to 287 378 finding a slope. 379 380 It also create's a frounde file. 288 381 """ 289 382 dx = 0.05 … … 306 399 pro_instance.outputdir + basename +".sww", 307 400 points, 308 pro_instance.outputdir + "fslope_depth_" + end,309 pro_instance.outputdir + "fslope_velocity_x_" + end,310 pro_instance.outputdir + "fslope_velocity_y_" + end,311 pro_instance.outputdir + "fslope_stage_" + end,312 pro_instance.outputdir + "fslope_froude_" + end,401 pro_instance.outputdir + slope_tag + "slope_depth_" + end, 402 pro_instance.outputdir + slope_tag + "slope_velocity_x_" + end, 403 pro_instance.outputdir + slope_tag + "slope_velocity_y_" + end, 404 pro_instance.outputdir + slope_tag + "slope_stage_" + end, 405 pro_instance.outputdir + slope_tag + "slope_froude_" + end, 313 406 time_thinning=1) 314 407 315 408 def find_min_slopes(times, slope_locations, slopes, 316 break_times,anuga_break_times, band_offset):317 bands = break_times2bands( break_times, band_offset)318 319 waves = {}409 anuga_break_times, band_offset): 410 bands = break_times2bands(anuga_break_times, band_offset) 411 412 waves = [] 320 413 for i,_ in enumerate(bands[0:-1]): 321 id = "wave " + str(i)322 414 max_q, max_q_times = get_min_in_band(bands[i], bands[i+1], 323 415 times, slopes) 324 waves[id] = {SLOPE_STR:max_q, 325 TIME_STR:max_q_times} 416 waves.append({SLOPE_STR:max_q, TIME_STR:max_q_times}) 326 417 return waves 327 418 … … 366 457 time_vector, quantity_array, 0) 367 458 #print "band_quantity",band_quantity 368 max_quantity_indices = argmin(band_quantity, axis=0) 459 try: 460 max_quantity_indices = argmin(band_quantity, axis=0) 461 except: 462 #print "time_vector", time_vector 463 print "min_time",min_time 464 print "max_time", max_time 465 return [],[] 466 369 467 #print "max_quantity_indices", max_quantity_indices 370 468 max_quantity_times = choose(max_quantity_indices, band_time) … … 386 484 bands.append(break_times[0]-0.5*(break_times[1]-break_times[0])) 387 485 486 388 487 for i,break_x in enumerate(break_times[0:-1]): 389 488 bands.append(break_times[i]+0.5*(break_times[i+1]-break_times[i])) … … 394 493 return bands 395 494 495 def plot_foude_slope_stage(wave_file, 496 break_time, 497 break_x, 498 save_as=None, 499 plot_title="", 500 is_interactive=False, 501 break_type="", 502 use_axis=None): 503 """ 504 """ 505 from pylab import ion, plot, xlabel, ylabel, close, legend, \ 506 savefig, title, axis, setp, subplot, grid, axvspan 507 from anuga.shallow_water.data_manager import csv2dict 508 509 510 511 # Load in the csv files and convert info from strings to floats 512 simulation, _ = csv2dict(wave_file) 513 location = [float(x) for x in simulation['x location']] 514 slope = [float(x) for x in simulation['min slope']] 515 time = [float(x) for x in simulation['Time']] 516 froude = [float(x) for x in simulation['Froude']] 517 518 min_location = min(location) 519 max_location = max(location) 520 521 if is_interactive: 522 ion() 523 # The upper subplot 524 subplot(311) 525 l_froude = plot(location, froude) 526 #setp(l_froude, color='r') 527 528 # Add axis stuff 529 title(plot_title) 530 y_label = "Froude Number" 531 ylabel(y_label) 532 grid(True) 533 axvspan(break_x-0.001,break_x+0.001, facecolor='g') 534 axis(ymin=0, ymax=1.8) 535 536 # The slope subplot 537 subplot(312) 538 l_slope = plot(location, slope) 539 setp(l_slope, color='r') 540 541 # Add axis stuff and legend 542 x_label = "X location, m" 543 y_label = "Stage slope" 544 #xlabel(x_label) 545 ylabel(y_label) 546 grid(True) 547 axvspan(break_x-0.001,break_x+0.001, facecolor='g') 548 axis(ymin=-0.5, ymax=0) 549 550 # The time, x location subplot 551 subplot(313) 552 l_time = plot(location, time) 553 setp(l_time, color='g') 554 #print "break_x", break_x 555 #print "break_time", break_time 556 plot([break_x], [break_time], 'yo') 557 #plot([break_x-1], [], 'yo') 558 559 # Add axis stuff and legend 560 x_label = "X location, m" 561 y_label = "time, sec" 562 xlabel(x_label) 563 ylabel(y_label) 564 grid(True) 565 566 567 # The order defines the label 568 #legend((legend_exp, legend_sim),'upper left') 569 #legend(('Wave front'),'upper left') 570 if use_axis is not None: 571 axis(use_axis) 572 573 if is_interactive: 574 # Wait for enter pressed 575 raw_input() 576 577 if save_as is not None: 578 savefig(save_as) 579 580 #Need to close this plot 581 close() 582 396 583 #------------------------------------------------------------- 397 584 if __name__ == "__main__": … … 400 587 from scenarios import scenarios 401 588 #scenarios = [scenarios[0]] 402 outputdir_tag = "_good_tri_area_0.01_limiterC" 589 outputdir_tag = "_good_tri_area_0.01_limiterD" 590 slope_tag = "" 591 #outputdir_tag = "_good_tri_area_0.01_limiterC" 592 #slope_tag = "f" 403 593 #outputdir_tag = "_test_limiterC" 404 594 #scenarios = [scenarios[0]] # !!!!!!!!!!!!!!!!!!!!!! 405 #gauges_for_slope(outputdir_tag, scenarios) 406 auto_graph_slopes(outputdir_tag, scenarios) #, is_interactive=True) 407 #auto_find_min_slopes(outputdir_tag, scenarios) 408 auto_graph_froudes(outputdir_tag, scenarios) 595 #scenarios = scenarios[4:] # !!!!!!!!!!!!!!!!!!!!!! 596 597 #gauges_for_slope(slope_tag, outputdir_tag, scenarios) 598 #auto_graph_slopes(outputdir_tag, scenarios) #, is_interactive=True) 599 #auto_find_min_slopes(slope_tag, outputdir_tag, scenarios) 600 #auto_graph_froudes(outputdir_tag, scenarios) 601 auto_plot_froude_slopes(slope_tag, outputdir_tag, scenarios) -
anuga_work/development/Hinwood_2008/test_slope.py
r5494 r5503 38 38 actual = array([0,0.01,0.02,0.03]) 39 39 self.assert_ (allclose(dtimes, actual, 0.001)) 40 actual = array([0.5, 0.50 1, 0.502])40 actual = array([0.5, 0.5005, 0.5015]) 41 41 self.assert_ (allclose(depths, actual, 0.001)) 42 42 … … 111 111 self.assert_ (allclose(ensure_numeric(bands), actual, 0.001)) 112 112 113 def test_get_time_band(self):113 def not_used_test_get_time_band(self): 114 114 time = array([0,1,2,3,4,5]) 115 115 slopes = array([[10,12], … … 158 158 [10,20,30,40]]) 159 159 self.assert_ (allclose(band_quantity, actual, 0.001)) 160 161 def test_find_froude(self): 162 times_froude = array([0,1]) 163 locations_froude = array([0,1,2,3,4,5]) 164 froudes = array([[0,1,2,3,4,5], 165 [0,10,20,30,40,50]]) 166 times = array([0,0,1,1]) 167 locations = array([0.5, 1.5, 2.5, 3.5]) 168 froudes = find_froude(times_froude, locations_froude, 169 froudes, times, locations) 170 actual = array([0.5, 1.5, 25, 35]) 171 self.assert_ (allclose(ensure_numeric(froudes), actual, 0.001)) 172 173 160 174 161 175 if __name__ == "__main__": 162 176 163 177 suite = unittest.makeSuite(slopeTestCase,'test') 164 suite = unittest.makeSuite(slopeTestCase,'test_get_band')178 #suite = unittest.makeSuite(slopeTestCase,'test_find_froude') 165 179 166 180 runner = unittest.TextTestRunner() #verbosity=2)
Note: See TracChangeset
for help on using the changeset viewer.
