Changeset 5503


Ignore:
Timestamp:
Jul 15, 2008, 5:00:20 PM (11 years ago)
Author:
duncan
Message:

Current Hinwood scenario - added plotting of froude number, stage slope and time vs location for each wave

Location:
anuga_work/development/Hinwood_2008
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/Hinwood_2008/run_dam.py

    r5494 r5503  
    6161         friction=0.012,  # planed wood. http://www.lmnoeng.com/manningn.htm
    6262         outputdir_name=None,
    63          run_type=0):
     63         run_type=0,
     64         end_tag = '_limiterD'):
    6465
    6566   
    6667    basename = 'zz_' + metadata_dic['scenario_id']
    67     end_tag = '_limiterD'
     68   
    6869    if run_type == 1:
    6970        outputdir_name += '_test' + end_tag
     
    261262    # 4 is 0.01
    262263    # 5 is 0.001
    263     run_type = 4
     264    run_type = 5
    264265    #for run_data in [scenarios[5]]:
    265266    #scenarios = scenarios[2:]
     
    269270                             run_data,
    270271                             run_type = run_type,
    271                              outputdir_name=run_data['scenario_id'])
     272                             outputdir_name=run_data['scenario_id'],
     273                             end_tag='_limiterD')
    272274        #gauges_for_slope(pro_instance.outputdir,[run_data])
  • anuga_work/development/Hinwood_2008/slope.py

    r5495 r5503  
    1010
    1111from Numeric import arange, array, zeros, Float, where, greater, less, \
    12      compress, argmin, choose
     12     compress, argmin, choose, searchsorted
    1313
    1414from anuga.fit_interpolate.interpolate import interpolate_sww2csv
     
    6969   
    7070    for i in arange(n_slope_locations):
    71         slope_locations[i] = (locations[i+1] + locations[i+1])/2.
    7271        delta_locations[i] = (locations[i+1] - locations[i])
     72        slope_locations[i] = locations[i] + 0.5*delta_locations[i]
    7373   
    7474    for j in arange(n_time):
     
    118118
    119119    # Add the contour line levels to the colorbar
    120     cbar.ax.set_ylabel('verbosity coefficient')
     120    cbar.ax.set_ylabel('stage slope')
    121121    #cbar.add_lines(CS2)
    122122         
     
    146146       
    147147    # 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,
    151155                  origin=origin)
    152156   
     
    168172
    169173    # Add the contour line levels to the colorbar
    170     cbar.ax.set_ylabel('verbosity coefficient')
     174    cbar.ax.set_ylabel('Froude Number')
    171175    #cbar.add_lines(CS2)
    172176         
     
    190194            anuga_break_times.append( \
    191195                break_time -  run_data['ANUGA_start_time'])
    192         plot_title = "Stage slope " + id
    193196        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
    194200        save_as = pro_instance.plots_dir + sep + \
    195201                            outputdir_name + "_slope_stage" + plot_type
    196202        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)
    202208        graph_contours(times, locations, slopes,
    203209                       plot_title=plot_title,
     
    221227            anuga_break_times.append( \
    222228                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
    225231        froude_file = pro_instance.outputdir + "fslope_froude_" + end
     232        print "Creating ", froude_file
    226233        save_as = pro_instance.plots_dir + sep + \
    227234                            outputdir_name + "_froude" + plot_type
     
    229236        dtimes, sensors = get_band(anuga_break_times[0]-TIME_BORDER,
    230237                                   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)
    233241        #print "dtimes", dtimes
    234242        #print "sensors", sensors
     
    241249                     is_interactive=is_interactive)
    242250       
    243  
    244 def auto_find_min_slopes(outputdir_tag, scenarios):
     251def 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   
     281def auto_find_min_slopes(slope_tag, outputdir_tag, scenarios):
    245282    """
    246283   
     
    258295                break_time -  run_data['ANUGA_start_time'])
    259296           
    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
    261300       
    262301        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
    263305        waves = find_min_slopes(times, slope_locations, slopes,
    264                                 run_data['break_times'],
    265306                                anuga_break_times,
    266307                                run_data['band_offset'])
     308       
    267309        # 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"
    271315            print "wave_file", wave_file
     316            froudes = find_froude(times_froude, locations_froude,
     317                             froudes_a, wave[TIME_STR],
     318                                  slope_locations)
    272319            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"])
    274321            wave_writer.writerows(map(None,
    275322                                      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
    278327           
    279             #wave_writer.close()
     328def 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'])
    280350           
     351        #run_data['break_type'] = (run_data['break_type'][0])
     352        for i in range(len(run_data['break_type'])):
    281353           
    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   
     375def gauges_for_slope(slope_tag, outputdir_tag, scenarios):
    285376    """
    286377    This is used to create a stage file, using gauges relivent to
    287378    finding a slope.
     379
     380    It also create's a frounde file.
    288381    """
    289382    dx = 0.05
     
    306399            pro_instance.outputdir + basename +".sww",
    307400            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,
    313406            time_thinning=1)
    314407
    315408def 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 = []
    320413    for i,_ in enumerate(bands[0:-1]):
    321         id = "wave " + str(i)
    322414        max_q, max_q_times = get_min_in_band(bands[i], bands[i+1],
    323415                                             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})
    326417    return waves
    327418
     
    366457                                         time_vector, quantity_array, 0)
    367458    #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       
    369467    #print "max_quantity_indices", max_quantity_indices
    370468    max_quantity_times = choose(max_quantity_indices, band_time)
     
    386484    bands.append(break_times[0]-0.5*(break_times[1]-break_times[0]))
    387485
     486
    388487    for i,break_x in enumerate(break_times[0:-1]):
    389488        bands.append(break_times[i]+0.5*(break_times[i+1]-break_times[i]))
     
    394493    return bands
    395494
     495def 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   
    396583#-------------------------------------------------------------
    397584if __name__ == "__main__":
     
    400587    from scenarios import scenarios
    401588    #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"
    403593    #outputdir_tag = "_test_limiterC"
    404594    #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  
    3838        actual = array([0,0.01,0.02,0.03])
    3939        self.assert_ (allclose(dtimes, actual, 0.001))
    40         actual = array([0.5, 0.501, 0.502])
     40        actual = array([0.5, 0.5005, 0.5015])
    4141        self.assert_ (allclose(depths, actual, 0.001))
    4242
     
    111111        self.assert_ (allclose(ensure_numeric(bands), actual, 0.001))
    112112       
    113     def test_get_time_band(self):
     113    def not_used_test_get_time_band(self):
    114114        time = array([0,1,2,3,4,5])
    115115        slopes = array([[10,12],
     
    158158                        [10,20,30,40]])
    159159        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       
    160174       
    161175if __name__ == "__main__":
    162176
    163177    suite = unittest.makeSuite(slopeTestCase,'test')
    164     suite = unittest.makeSuite(slopeTestCase,'test_get_band')
     178    #suite = unittest.makeSuite(slopeTestCase,'test_find_froude')
    165179
    166180    runner = unittest.TextTestRunner()  #verbosity=2)
Note: See TracChangeset for help on using the changeset viewer.