"""Read in sww file, interpolate at specified locations and Cross section at specified T """ from os import sep import Numeric import project from pyvolution.util import file_function from create_buildings import WidtH # from run_friction import fric from pylab import * fric = 'mid_limit_0' # FILE LABEL #swwfile = project.newoutputname + '.sww' #swwfile = project.outputname #swwfile_B = project.outputdir + sep + 'Buildings_3790.sww' swwfile_B = project.outputdir + sep + 'buildings_Str_D=5.8_976.sww' #swwfile_F = project.outputdir + sep + 'buildings_Str_half_Off_D=6.89_1116.sww' #swwfile_B = project.outputdir + sep + 'friction_n=6_3135.sww' swwfile_F = project.outputdir + sep + 'friction_n=0.9_3135.sww' # gauge_depth = Numeric.arrayrange(51, 1.5*WidtH, 25) # Random special #gauge_depth = Numeric.arrayrange(0, 5*WidtH, 5) gauge_depth = [49.600000000000001, 74.6, 99.6, 124.6, 149.6, 174.6, 199.6, 224.6, 240.00999999999999] gauge_breadth = 100 gauge_locations = [] for GD in gauge_depth: gauge_location = [GD,gauge_breadth] gauge_locations.append(gauge_location) #Read model output quantities = ['stage', 'xmomentum', 'ymomentum'] f_B = file_function(swwfile_B, quantities = quantities, interpolation_points = gauge_locations, verbose = True, use_cache = False) f_F = file_function(swwfile_F, quantities = quantities, interpolation_points = gauge_locations, verbose = True, use_cache = False) T = Numeric.arrayrange(700,1000,10) #T = [ 20, 50, 100, 150, 200 ] for t in T: model_time = [] stages_B = [] stages_F = [] elevations = [] momenta = [] velocity = [] for i, GL in enumerate(gauge_depth): # mines data from specified time values in sww file wB = f_B(t, point_id = i)[0] uhB = f_B(t, point_id = i)[1] vhB = f_B(t, point_id = i)[2] wF = f_F(t, point_id = i)[0] uhF = f_B(t, point_id = i)[1] vhF = f_B(t, point_id = i)[2] #m = sqrt(uh*uh + vh*vh) #Absolute momentum velB = sqrt(uhB*uhB + vhB*vhB) / (wB + 1.e-30) #Absolute velocity velF = sqrt(uhF*uhF + vhF*vhF) / (wF + 1.e-30) #Absolute velocity #print velB, "buiding velocity" #print velF, "friction velocity" stages_B.append(wB) stages_F.append(wF) #elevations.append(zB) #Should be constant over time #momenta.append(mB) #velocityB.append(velB) #velocityF.append(velF) #print len(stages), "stages" diff = abs(Numeric.array(stages_B) - Numeric.array(stages_F)) #print diff, "difference" ion() hold(True) subplot(211) plot(gauge_depth, stages_B, '-b',gauge_depth,stages_F,'-r') title('Stage@Time_%d_ vs Gauge length (m)[red: fric, blue: buildings] ' %t) xlabel('gauge length(m)') ylabel('water depths (m)') #savefig('Comp_Time%d_B_vs_F' %t) hold(False) print stages_F, " stages" print gauge_depth, " gauge depths" diff_area_F = trapz(gauge_depth,stages_F) diff_area_B = trapz(gauge_depth,stages_B) diff_area = abs(diff_area_F) - abs(diff_area_B) print diff_area, " <<< diff area*********" #savefig('Final_comp_n=%s' %fric) subplot(212) subplots_adjust(hspace = 0.4) plot(gauge_depth, diff, 'g') title('Final Differences in water levels') xlabel('gauge length(m)') ylabel('water differences (m)') savefig('Final_Water_diff_n=%s' %fric) print "finished"