1 | """Read in sww file, interpolate at specified locations and |
---|
2 | Cross section at specified T. Cross sections are spatial drawing |
---|
3 | at the specified time. |
---|
4 | |
---|
5 | '2nd normal' test is employed to find the best fit friction value |
---|
6 | for the given building scenario => 'B_file' |
---|
7 | |
---|
8 | """ |
---|
9 | from os import sep |
---|
10 | import Numeric |
---|
11 | import project |
---|
12 | from pyvolution.util import file_function |
---|
13 | from create_buildings import WidtH |
---|
14 | from pylab import * |
---|
15 | from numerical_tools import norm, corr, err |
---|
16 | |
---|
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'] |
---|
25 | |
---|
26 | |
---|
27 | fd = open("buildings_Orth(45)_report.csv", 'w') |
---|
28 | fd.write( 'buildings_Orth(45)_report. Summary at bottom.' + '\n' ) |
---|
29 | |
---|
30 | |
---|
31 | # Building scenario file to be compared to |
---|
32 | #--------------------------------------------- |
---|
33 | #B_list = ['buildings_S0_D=23_1784.sww'] |
---|
34 | #--------------------------------------------- |
---|
35 | |
---|
36 | best_fric=[] |
---|
37 | B_Density=[] |
---|
38 | best_norm=[] |
---|
39 | Breadth=[] |
---|
40 | |
---|
41 | |
---|
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]) |
---|
47 | |
---|
48 | gauge_depth = Numeric.arrayrange((52.5-bre/2), 1.5*WidtH, 25) |
---|
49 | gauge_breadth = 100 |
---|
50 | gauge_locations = [] |
---|
51 | |
---|
52 | for GD in gauge_depth: |
---|
53 | gauge_location = [GD,gauge_breadth] |
---|
54 | gauge_locations.append(gauge_location) |
---|
55 | |
---|
56 | |
---|
57 | #Read model output |
---|
58 | quantities = ['stage', 'xmomentum', 'ymomentum'] |
---|
59 | |
---|
60 | f_B = file_function(swwfile_B, |
---|
61 | quantities = quantities, |
---|
62 | interpolation_points = gauge_locations, |
---|
63 | verbose = True, |
---|
64 | use_cache = True) |
---|
65 | |
---|
66 | fric_OK=[] |
---|
67 | least_norm=[100] |
---|
68 | Normal=[] |
---|
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) |
---|
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) |
---|
80 | |
---|
81 | t=990 # specifies times to look at cross sections |
---|
82 | |
---|
83 | stages_B = [] |
---|
84 | stages_F = [] |
---|
85 | vel_B=[] |
---|
86 | vel_F=[] |
---|
87 | for i, GL in enumerate(gauge_depth): |
---|
88 | |
---|
89 | k = len(gauge_depth) # point of velocity comparision |
---|
90 | # mines data from specified time values in sww file |
---|
91 | # Building file |
---|
92 | wB = f_B(t, point_id = i)[0] |
---|
93 | uhB = f_B(t, point_id = i)[1] |
---|
94 | vhB = f_B(t, point_id = i)[2] |
---|
95 | # Friction file |
---|
96 | wF = f_F(t, point_id = i)[0] |
---|
97 | uhF = f_B(t, point_id = i)[1] |
---|
98 | vhF = f_B(t, point_id = i)[2] |
---|
99 | |
---|
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) |
---|
105 | |
---|
106 | |
---|
107 | vel_B.append(velB) |
---|
108 | vel_F.append(velF) |
---|
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 |
---|
120 | |
---|
121 | else: |
---|
122 | print fric, " NOT OK" |
---|
123 | |
---|
124 | |
---|
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) |
---|
131 | |
---|
132 | print vel_F, "friction" |
---|
133 | print vel_B, "Buildings" |
---|
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" |
---|
143 | |
---|
144 | |
---|
145 | len_file=Numeric.arrayrange(0,len(B_list),1) |
---|
146 | fd.write('_______________________________________________________' + '\n' ) |
---|
147 | fd.write( '<<<<< buildings_Orth(45)_report Summary >>>>>' + '\n' ) |
---|
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 | |
---|