[2417] | 1 | """Read in sww file, interpolate at specified locations and |
---|
[2458] | 2 | Cross section at specified T. Cross sections are spatial drawing |
---|
| 3 | at the specified time. |
---|
[2417] | 4 | |
---|
[2458] | 5 | '2nd normal' test is employed to find the best fit friction value |
---|
| 6 | for the given building scenario => 'B_file' |
---|
| 7 | |
---|
[2417] | 8 | """ |
---|
| 9 | from os import sep |
---|
| 10 | import Numeric |
---|
| 11 | import project |
---|
| 12 | from pyvolution.util import file_function |
---|
[2458] | 13 | from create_buildings import WidtH |
---|
[2417] | 14 | from pylab import * |
---|
[2429] | 15 | from numerical_tools import norm, corr, err |
---|
[2417] | 16 | |
---|
[2466] | 17 | B_list=['buildings_Orth(45)_D=3_1566.sww', |
---|
| 18 | 'buildings_Orth(45)_D=5_1247.sww', |
---|
| 19 | 'buildings_Orth(45)_D=7_1270.sww', |
---|
| 20 | 'buildings_Orth(45)_D=9_804.sww', |
---|
| 21 | 'buildings_Orth(45)_D=11_844.sww', |
---|
| 22 | 'buildings_Orth(45)_D=13_918.sww', |
---|
| 23 | 'buildings_Orth(45)_D=15_1250.sww', |
---|
| 24 | 'buildings_Orth(45)_D=17_1735.sww'] |
---|
[2417] | 25 | |
---|
[2418] | 26 | |
---|
[2466] | 27 | fd = open("buildings_Orth(45)_report.csv", 'w') |
---|
| 28 | fd.write( 'buildings_Orth(45)_report. Summary at bottom.' + '\n' ) |
---|
[2446] | 29 | |
---|
[2417] | 30 | |
---|
[2458] | 31 | # Building scenario file to be compared to |
---|
| 32 | #--------------------------------------------- |
---|
| 33 | #B_list = ['buildings_S0_D=23_1784.sww'] |
---|
| 34 | #--------------------------------------------- |
---|
[2417] | 35 | |
---|
[2458] | 36 | best_fric=[] |
---|
| 37 | B_Density=[] |
---|
| 38 | best_norm=[] |
---|
| 39 | Breadth=[] |
---|
[2417] | 40 | |
---|
[2446] | 41 | |
---|
[2458] | 42 | for B_file in B_list: |
---|
| 43 | swwfile_B = project.outputdir + sep + B_file |
---|
| 44 | a=B_file.split('_') |
---|
| 45 | b=a[2].split('=') |
---|
| 46 | bre=int(b[1]) |
---|
[2417] | 47 | |
---|
[2458] | 48 | gauge_depth = Numeric.arrayrange((52.5-bre/2), 1.5*WidtH, 25) |
---|
| 49 | gauge_breadth = 100 |
---|
| 50 | gauge_locations = [] |
---|
[2429] | 51 | |
---|
[2458] | 52 | for GD in gauge_depth: |
---|
| 53 | gauge_location = [GD,gauge_breadth] |
---|
| 54 | gauge_locations.append(gauge_location) |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | #Read model output |
---|
[2466] | 58 | quantities = ['stage', 'xmomentum', 'ymomentum'] |
---|
[2458] | 59 | |
---|
| 60 | f_B = file_function(swwfile_B, |
---|
[2418] | 61 | quantities = quantities, |
---|
| 62 | interpolation_points = gauge_locations, |
---|
| 63 | verbose = True, |
---|
| 64 | use_cache = True) |
---|
[2417] | 65 | |
---|
[2458] | 66 | fric_OK=[] |
---|
| 67 | least_norm=[100] |
---|
| 68 | Normal=[] |
---|
[2466] | 69 | Friction_large = list(Numeric.arange(1,120,1)) |
---|
| 70 | Friction = list(Numeric.arange(0.025,1,0.025)) |
---|
| 71 | #Friction_log = [0.025, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000] |
---|
| 72 | Friction.extend(Friction_large) |
---|
[2458] | 73 | for fric in Friction: |
---|
| 74 | swwfile_F = project.outputdir + sep + 'friction_n=' + str(fric) + '_3135.sww' |
---|
| 75 | f_F = file_function(swwfile_F, |
---|
| 76 | quantities = quantities, |
---|
| 77 | interpolation_points = gauge_locations, |
---|
| 78 | verbose = True, |
---|
| 79 | use_cache = True) |
---|
[2417] | 80 | |
---|
[2458] | 81 | t=990 # specifies times to look at cross sections |
---|
[2417] | 82 | |
---|
[2458] | 83 | stages_B = [] |
---|
| 84 | stages_F = [] |
---|
[2466] | 85 | vel_B=[] |
---|
| 86 | vel_F=[] |
---|
[2458] | 87 | for i, GL in enumerate(gauge_depth): |
---|
[2417] | 88 | |
---|
[2466] | 89 | k = len(gauge_depth) # point of velocity comparision |
---|
[2458] | 90 | # mines data from specified time values in sww file |
---|
| 91 | # Building file |
---|
| 92 | wB = f_B(t, point_id = i)[0] |
---|
[2466] | 93 | uhB = f_B(t, point_id = i)[1] |
---|
| 94 | vhB = f_B(t, point_id = i)[2] |
---|
[2458] | 95 | # Friction file |
---|
| 96 | wF = f_F(t, point_id = i)[0] |
---|
[2466] | 97 | uhF = f_B(t, point_id = i)[1] |
---|
| 98 | vhF = f_B(t, point_id = i)[2] |
---|
[2417] | 99 | |
---|
[2458] | 100 | #m = sqrt(uh*uh + vh*vh) #Absolute momentum |
---|
| 101 | velB = sqrt(uhB*uhB + vhB*vhB) / (wB + 1.e-30) #Absolute velocity |
---|
| 102 | velF = sqrt(uhF*uhF + vhF*vhF) / (wF + 1.e-30) #Absolute velocity |
---|
| 103 | stages_B.append(wB) |
---|
| 104 | stages_F.append(wF) |
---|
[2466] | 105 | |
---|
| 106 | |
---|
| 107 | vel_B.append(velB) |
---|
| 108 | vel_F.append(velF) |
---|
[2458] | 109 | Ar_B = Numeric.array(stages_B) |
---|
| 110 | Ar_F = Numeric.array(stages_F) |
---|
| 111 | diff = abs(Ar_B - Ar_F) |
---|
| 112 | norm_F = err(Ar_F, Ar_B, 2, relative = True) # 2nd norm (rel. RMS) test |
---|
| 113 | Normal.append(norm_F) |
---|
| 114 | #print " " |
---|
| 115 | #print norm_F, "2nd norm test result (N_F)" |
---|
| 116 | #print " " |
---|
| 117 | if norm_F < least_norm: |
---|
| 118 | least_norm = norm_F |
---|
| 119 | fric_OK = fric |
---|
[2466] | 120 | |
---|
[2458] | 121 | else: |
---|
| 122 | print fric, " NOT OK" |
---|
[2466] | 123 | |
---|
[2417] | 124 | |
---|
[2458] | 125 | print fric_OK, " <- This Manning's n is best fit from 2nd norm test" |
---|
| 126 | print least_norm, "2nd normal of Friction vs Buildings " |
---|
| 127 | best_fric.append(fric_OK) |
---|
| 128 | B_Density.append((float(bre)**2)/625) |
---|
| 129 | best_norm.append(least_norm) |
---|
| 130 | Breadth.append(bre) |
---|
[2429] | 131 | |
---|
[2466] | 132 | print vel_F, "friction" |
---|
| 133 | print vel_B, "Buildings" |
---|
[2458] | 134 | fd.write('_______________________________________________________' + '\n' ) |
---|
| 135 | fd.write( '<<<< Mannings n for [' + B_file + '] >>>>' + '\n' ) |
---|
| 136 | fd.write( '(n), 2nd norm ' + '\n') |
---|
| 137 | for j,drt in enumerate(Friction): |
---|
| 138 | fd.write(str(Friction[j]) + ', ' + str(Normal[j]) + '\n') |
---|
| 139 | |
---|
| 140 | print best_fric, " << Best fit friction" |
---|
| 141 | print B_Density, " << Building footprint density" |
---|
| 142 | print best_norm, " << Matching normal densities" |
---|
[2417] | 143 | |
---|
[2458] | 144 | |
---|
| 145 | len_file=Numeric.arrayrange(0,len(B_list),1) |
---|
| 146 | fd.write('_______________________________________________________' + '\n' ) |
---|
[2466] | 147 | fd.write( '<<<<< buildings_Orth(45)_report Summary >>>>>' + '\n' ) |
---|
[2458] | 148 | fd.write( 'width, footprint ratio, mannings n, 2nd norm ' + '\n') |
---|
| 149 | for i in len_file: |
---|
| 150 | fd.write( str(Breadth[i]) + ', ' + str(B_Density[i]) + ', ' + str(best_fric[i]) + ', ' + str(best_norm[i]) + '\n') |
---|
| 151 | fd.close |
---|
| 152 | |
---|