source: development/momentum_sink/CCS.py @ 2465

"""Read in sww file, interpolate at specified locations and
    Cross section at specified T. Cross sections are spatial drawing
    at the specified time.

    '2nd normal' test is employed to find the best fit friction value
    for the given building scenario => 'B_file'

"""
from os import sep
import Numeric
import project
from pyvolution.util import file_function
from create_buildings import WidtH
from pylab import *
from numerical_tools import norm, corr, err  

B_list=['buildings_Str_D=3_1714.sww',
        'buildings_Str_D=5_1578.sww',
        'buildings_Str_D=7_1067.sww',
        'buildings_Str_D=9_1052.sww',
        'buildings_Str_D=11_694.sww',
        'buildings_Str_D=13_693.sww',
        'buildings_Str_D=15_717.sww',
        'buildings_Str_D=17_749.sww',
        'buildings_Str_D=19_1008.sww',
        'buildings_Str_D=21_1581.sww',
        'buildings_Str_D=23_1798.sww']


fd = open("Straight_log_report.csv", 'w')
fd.write( 'Straight_log_report. Summary at bottom.' + '\n' )


# Building scenario file to be compared to
#---------------------------------------------
#B_list = ['buildings_S0_D=23_1784.sww']
#---------------------------------------------

best_fric=[]
B_Density=[]
best_norm=[]
Breadth=[]


for B_file in B_list:
    swwfile_B = project.outputdir + sep  + B_file
    a=B_file.split('_')
    b=a[2].split('=')
    bre=int(b[1])

    gauge_depth = Numeric.arrayrange((52.5-bre/2), 1.5*WidtH, 25)
    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', 'elevation', 'xmomentum', 'ymomentum']

    f_B = file_function(swwfile_B,
                      quantities = quantities,
                      interpolation_points = gauge_locations,
                      verbose = True,
                      use_cache = True)

    fric_OK=[]
    least_norm=[100]
    Normal=[]
    #Friction = Numeric.arrayrange(1,120,1)
    Friction = [0.025, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]
    for fric in Friction:
        swwfile_F = project.outputdir + sep  + 'friction_n=' + str(fric) + '_3135.sww'
        f_F = file_function(swwfile_F,
                          quantities = quantities,
                          interpolation_points = gauge_locations,
                          verbose = True,
                          use_cache = True)

        t=990 # specifies times to look at cross sections
        
        stages_B = []
        stages_F = []
        for i, GL in enumerate(gauge_depth):

            # mines data from specified time values in sww file
            # Building file
            wB = f_B(t, point_id = i)[0]
            zB = f_B(t, point_id = i)[1]
            uhB = f_B(t, point_id = i)[2]
            vhB = f_B(t, point_id = i)[3]
            # Friction file
            wF = f_F(t, point_id = i)[0]        
            uhF = f_B(t, point_id = i)[2]
            vhF = f_B(t, point_id = i)[3]

            #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                   
            stages_B.append(wB)
            stages_F.append(wF)
            
            
        Ar_B = Numeric.array(stages_B)
        Ar_F = Numeric.array(stages_F)
        diff = abs(Ar_B - Ar_F)
        norm_F = err(Ar_F, Ar_B, 2, relative = True)  # 2nd norm (rel. RMS) test
        Normal.append(norm_F)
        #print " "
        #print norm_F, "2nd norm test result (N_F)"
        #print " "
        if norm_F < least_norm:
            least_norm = norm_F
            fric_OK = fric
        else:
            print fric, " NOT OK" 
    
    print fric_OK, " <- This Manning's n is best fit from 2nd norm test"
    print least_norm, "2nd normal of Friction vs Buildings  "
    best_fric.append(fric_OK)
    B_Density.append((float(bre)**2)/625)
    best_norm.append(least_norm)
    Breadth.append(bre)
    
    print len(Friction), "friction"
    print len(Normal), "normal"
    fd.write('_______________________________________________________' + '\n' )
    fd.write( '<<<< Mannings n for [' + B_file + '] >>>>' + '\n' )
    fd.write( '(n), 2nd norm ' + '\n')
    for j,drt in enumerate(Friction):
        fd.write(str(Friction[j]) + ', ' +  str(Normal[j]) + '\n')
         
print best_fric, " << Best fit friction"
print B_Density, " << Building footprint density"
print best_norm, " << Matching normal densities"


len_file=Numeric.arrayrange(0,len(B_list),1)
fd.write('_______________________________________________________' + '\n' )
fd.write( '<<<<< Straight_log_report Summary >>>>>' + '\n' )
fd.write( 'width, footprint ratio, mannings n, 2nd norm ' + '\n')
for i in len_file:
    fd.write( str(Breadth[i]) + ', ' + str(B_Density[i]) + ', ' +  str(best_fric[i]) + ', ' +  str(best_norm[i]) + '\n')
fd.close