source: anuga_work/development/momentum_sink/scripts/CCS.py @ 6939

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

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

    the RMS result for each friction file is written to a csv file
    

"""
from os import sep
import os
import Numeric
import project
import string
from anuga.pyvolution.util import file_function
from numerical_tools import norm, corr, err  

# NAME FOLDER TO GET BUILDING SCENARIOS FROM
#-------------------------------------------------------------------------
folder_name = 'buildings_Rot(45)_10'   # <<<<<< place folder name here <<<<<<<
#-------------------------------------------------------------------------
wave_height = str(10)                    # <<<<<< place wave depth here <<<<<<<
#-------------------------------------------------------------------------



# FILE SORTER
#-------------------------------------------------------------------------
# Sorting function used to sort files in the specified folder #    
B_list_unsorted=os.listdir(project.building_dir+ folder_name + sep) 
B_strings=[]
B_list=[]
for B_L in B_list_unsorted:
    A = B_L.split('=')[0] + '=' # first part of string
    B = int(B_L.split('=')[1].split('_')[0]) # middle part of string
    C = '_' + B_L.split('=')[1].split('_')[1] # last part of string
    B_current = [A,B,C]; B_strings.append(B_current)

B_strings.sort(lambda x, y: cmp(x[1],y[1]))

for B_str in B_strings:
    B_str[1]=str(B_str[1])
    final=string.join(B_str).replace(' ','')
    B_list.append(final)
print B_list
print "-----------------------------------------------------"
print ">> ",B_list[0],"<< FIRST file"
print ">> ",B_list[-1],"<< LAST file"
print "-----------------------------------------------------"
# End of sorting file

#-------------------------------------------------------------------------
print "-----------------------------------------------------"
fd = open(folder_name + '_report.csv', 'w')
fd.write( folder_name + '_report from CCS.py using 2nd normal test' + '\n' )

#-------------------------------------------------------------------------
# Specify correct building parameters of scenario. Only used in naming 
WidtH = 200
BL = 25

fd.write( 'Building scenario parameters: WidtH = ' + str(WidtH) + ', Block length = ' + str(BL) + ', Wave height = ' + str(wave_height) + '\n' )
print " "
print "-----------------------------------------------------"
print 'Building scenario parameters: WidtH = ' + str(WidtH) + ', Block length = ' + str(BL) + ', Wave height = ' + str(wave_height)
print "-----------------------------------------------------"
#-------------------------------------------------------------------------

best_fric=[]
B_Density=[]
best_norm=[]
Breadth=[]
prev_index=0
u=1

#-------------------------------------------------------------------------
# Buildings scenario loops from sorting. loops through sww files in specifies folder name
for B_file in B_list:
    swwfile_B = project.building_dir + folder_name +sep+ B_file
    bre = float(B_file.split('=')[1].split('_')[0]) # middle part of string
    print bre, "breadth"
    
    # IMPORTANT, locates gauge points between buildings.
    gauge_depth = list(Numeric.arange(0.2*WidtH + ((BL-bre)/2), 1.2*WidtH+0.1, 25))
    gauge_depth.append((1.2*WidtH)+0.01)
    gauge_breadth = 100 # identifies where the cross section will be cut
    print gauge_depth, "gauge locations"
    # raw_input('Check gauge locations >> are they in order?')
    
    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)

    fd.write('_______________________________________________________' + '\n' )
    fd.write( '<<<< Mannings n for [' + B_file + '] >>>>' + '\n' )
    fd.write( '(n), 2nd norm ' + '\n')

    fric_OK=[]
    least_norm=[100]
    Normal=[]

    # creation of list of friction values to check through, make sure sww files exist before proceeding.
    Friction_large = list(Numeric.arange(1,120,1))
    Friction = list(Numeric.arange(0.025,1,0.025))
    Friction_log = list([200, 500, 1000, 2000, 5000, 10000])
    Friction.extend(Friction_large)
    Friction.extend(Friction_log)
    #Friction_log = [0.025, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 1000, 2000, 5000, 10000]
    
    if prev_index < 6:
        pri=0
    else:
        pri=prev_index-5
        
    #-------------------------------------------------------------------------
   

    for f_index, fric in enumerate(Friction):
    #for f_index, fric in enumerate(Friction[pri:]):
        swwfile_F = project.friction_dir +'friction_files_'+wave_height+ sep + 'friction_' +wave_height+ '_n=' + str(fric) + '_3125.sww'
        f_F = file_function(swwfile_F,
                          quantities = quantities,
                          interpolation_points = gauge_locations,
                          verbose = True,
                          use_cache = False)
        
        t=990 # specifies times to look at cross sections
        stages_B = []
        stages_F = []
        vel_B=[]
        vel_F=[]
        for i, GL in enumerate(gauge_depth):

            # mines data from specified time values in sww file

            # Building values
            wB = f_B(t, point_id = i)[0]
            uhB = f_B(t, point_id = i)[1]
            vhB = f_B(t, point_id = i)[2]
            # Friction values
            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 = ((uhB*uhB + vhB*vhB)**0.5) / (wB + 1.e-30) #Absolute velocity
            velF = ((uhF*uhF + vhF*vhF)**0.5) / (wF + 1.e-30) #Absolute velocity                   
            stages_B.append(wB)
            stages_F.append(wF)

        vel_B.append(velB)
        vel_F.append(velF)
        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)
        fd.write(str(fric) + ', ' +  str(norm_F) + '\n')

        # finds lowest RMS thus friction value 
        if norm_F < least_norm:
            least_norm = norm_F
            fric_OK = fric
            prev_index = Friction.index(fric)
            
        else:
            print "not OK"
            #break

    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 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( '<<<<< ' + folder_name + ' 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