source: development/momentum_sink/CCS.py @ 2466

Last change on this file since 2466 was 2466, checked in by nicholas, 18 years ago

Full set of friction values run through CCS.py for various building layouts. Report started for produced reports in Experiment 1.doc

File size: 5.1 KB
Line 
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"""
9from os import sep
10import Numeric
11import project
12from pyvolution.util import file_function
13from create_buildings import WidtH
14from pylab import *
15from numerical_tools import norm, corr, err 
16
17B_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
27fd = open("buildings_Orth(45)_report.csv", 'w')
28fd.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
36best_fric=[]
37B_Density=[]
38best_norm=[]
39Breadth=[]
40
41
42for 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         
140print best_fric, " << Best fit friction"
141print B_Density, " << Building footprint density"
142print best_norm, " << Matching normal densities"
143
144
145len_file=Numeric.arrayrange(0,len(B_list),1)
146fd.write('_______________________________________________________' + '\n' )
147fd.write( '<<<<< buildings_Orth(45)_report Summary >>>>>' + '\n' )
148fd.write( 'width, footprint ratio, mannings n, 2nd norm ' + '\n')
149for i in len_file:
150    fd.write( str(Breadth[i]) + ', ' + str(B_Density[i]) + ', ' +  str(best_fric[i]) + ', ' +  str(best_norm[i]) + '\n')
151fd.close
152
Note: See TracBrowser for help on using the repository browser.