Changeset 2458

Timestamp:

Mar 1, 2006, 1:59:25 PM (18 years ago)

Author:

nicholas

Message:

CCS.py modified to loop through building scenarios and write results to file.

Location:

development/momentum_sink

Files:

• CCS.py (modified) (3 diffs)
• PCS.py (modified) (3 diffs)
• create_buildings.py (modified) (6 diffs)
• friction_block.py (modified) (4 diffs)
• loop_buildings.py (modified) (5 diffs)
• loop_friction.py (modified) (2 diffs)
• run_buildings.py (modified) (4 diffs)

development/momentum_sink/CCS.py

@@ -1 +1 @@
 """Read in sww file, interpolate at specified locations and
-    Cross section at specified T
+    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'
 
 """
@@ -7 +11 @@
 import project
 from pyvolution.util import file_function
-from create_buildings import porosity, WidtH
-# from run_friction import fric
+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']
 
 
-
-gauge_depth = Numeric.arrayrange(46, 1.5*WidtH, 25)
-print len(gauge_depth),"gauge_depth length"
-gauge_breadth = 100
-gauge_locations = []
-
-for GD in gauge_depth:
-    gauge_location = [GD,gauge_breadth]
-    gauge_locations.append(gauge_location)
+fd = open("Straight_log_report.csv", 'w')
+fd.write( 'Straight_log_report. Summary at bottom.' + '\n' )
 
 
-swwfile_B = project.outputdir + sep  + 'Straight_offset_11.sww'
+# Building scenario file to be compared to
+#---------------------------------------------
+#B_list = ['buildings_S0_D=23_1784.sww']
+#---------------------------------------------
 
-#Read model output
-quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
+best_fric=[]
+B_Density=[]
+best_norm=[]
+Breadth=[]
 
-f_B = file_function(swwfile_B,
-                  quantities = quantities,
-                  interpolation_points = gauge_locations,
-                  verbose = True,
-                  use_cache = True)
 
-Corr_OK=[]
-fric_OK=[]
-least_norm=[1]
-norm_OK=[]
-Friction = Numeric.arrayrange(1,80,1)
-for fric in Friction:
-    swwfile_F = project.outputdir + sep  + 'friction_n=' + repr(fric) + '_3135.sww'
-    f_F = file_function(swwfile_F,
+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,
@@ -49 +67 @@
                       use_cache = True)
 
-    t=1000 # specifies times to look at cross sections
+    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" 
     
-    stages_B = []
-    stages_F = []
-    for i, GL in enumerate(gauge_depth):
+    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"
 
-        # 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_B, Ar_F, 2, relative = True)# 2nd norm (RMS) test
-     
-    print " "
-    print norm_F, "norm N_F"
-    print " "
-    norm_OK.append(norm_F)
+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
 
-    if norm_F < least_norm:
-        least_norm = norm_F
-        fric_OK = fric
-    else:
-        print fric, " NOT OK" 
-
-    
-    
-print " "
-print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
-print fric_OK, " <- This Manning's n is best fit from 2nd norm test"
-print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
-print least_norm, "2nd normal of Friction vs Buildings  "
-print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
-

development/momentum_sink/PCS.py

@@ -7 +7 @@
 import project
 from pyvolution.util import file_function
-from create_buildings import porosity, WidtH
+from create_buildings import  WidtH
 # from run_friction import fric
 from pylab import *
@@ -16 +16 @@
 
 #swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
-swwfile_B = project.outputdir + sep  + 'Straight_offset.sww'
+swwfile_B = project.outputdir + sep  + 'buildings_S0_D=19_1037.sww'
 #swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
-swwfile_F = project.outputdir + sep  + 'friction_n=51_3135.sww'
+swwfile_F = project.outputdir + sep  + 'friction_n=37_3135.sww'
 
-gauge_depth = Numeric.arrayrange(46, 1.5*WidtH, 25) # Random special
+gauge_depth = Numeric.arrayrange(51, 1.5*WidtH, 25) # Random special
 #gauge_depth = Numeric.arrayrange(0, 1.5*WidtH, 25)
 gauge_breadth = 100
@@ -90 +90 @@
     
     plot(gauge_depth, stages_B, '-b',gauge_depth,stages_F,'-r')
-    title('Time_%d_ vs water depth (n=%f, BP=%f)' %(t, fric, porosity))
+    title('Time_%d_ vs water depth ' %t)
     xlabel('gauge length(m)')
     ylabel('water depths (m)')

development/momentum_sink/create_buildings.py

@@ -7 +7 @@
 from coordinate_transforms.geo_reference import Geo_reference
 
-depth = 11 # depth of building side to oncoming wave
-breadth = 11 # breadth of building, width of building to oncoming wave
-#from loop_buildings import depth
-
 WidtH = 200 # width of boudary in metres
-#breadth = depth
-
-print "building footprint"
-print depth * breadth , "m^2"
-block = 625
-BL = block**0.5
-
-porosity = breadth/BL                           
-print porosity, " Building porosity"
 
 #random.uniform(0,1) #proper implemantation of random generator
 
-def create_mesh(maximum_triangle_area,
+def create_mesh(maximum_triangle_area, depth,
                 mesh_file=None,
                 triangles_in_name = False):
@@ -32 +19 @@
     """
 
-    
+    WidtH = 200 # width of boudary in metres
+    breadth = depth
+
+    print "building footprint"
+    print depth * breadth , "m^2"
+    block = 625
+    BL = block**0.5
+
+    porosity = breadth/BL                               
+    print porosity, " Building porosity"
     # create a mesh instance of class Mesh
     m = Mesh()
@@ -44 +40 @@
     whs = depth/2 
     lhs = breadth/2 
-    #Th = (45 *(3.142/180)) # sets an initial rotation      
-    Th = 0
+    Th = (45 *(3.142/180)) # sets an initial rotation      
+    #Th = 0
     ForDep = (0.2*WidtH) + (BL-whs)
     RearDep = 1.2*WidtH 
         
     Breadths = Numeric.arrayrange( -(BL/2), WidtH+(BL/2), (BL))
+    #Breadths = Numeric.arrayrange( (BL/2), WidtH-(BL/2), (BL)) # For ortho-offset
     print Breadths, "Breadths"
     Depths = Numeric.arrayrange( ForDep, RearDep, BL )
@@ -55 +52 @@
 
     for i,D in enumerate(Depths):
-        Breadths = Breadths + ((-1)**i)*(BL/2) #Used to offset buildings
+        #Breadths = Breadths + ((-1)**i)*(BL/2) #Used to offset buildings
         for B in Breadths:
             wh1 = (-whs) * math.cos(Th) + (-lhs) * math.sin(Th)
@@ -76 +73 @@
     else:
         if triangles_in_name is True:
-            mesh_file = mesh_file[:-4] + '_Ran_' + str(depth) + '_' + str(triangle_count) \
+            mesh_file = mesh_file[:-4] + '_Orth(45)_D=' + str(depth) + '_' + str(triangle_count) \
                         + mesh_file[-4:]
         m.export_mesh_file(mesh_file)
@@ -83 +80 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    _, triangle_count = create_mesh(10,mesh_file="test.tsh")
+    _, triangle_count = create_mesh(100,15,mesh_file="test.tsh")
     print "triangle_count",triangle_count 

development/momentum_sink/friction_block.py

@@ -6 +6 @@
 from coordinate_transforms.geo_reference import Geo_reference
 
-#test = 300
+
 def create_mesh(maximum_triangle_area,
                 mesh_file=None,
@@ -17 +17 @@
     m = Mesh()
     
-    #print "test", test  
     # Boundary of problem
     WidtH = 200 # width of boudary in metres
-    #W = WidtH/8
-    #L = W
-    #outer_polygon = [[0,0],[2*WidtH,0],[2*WidtH,WidtH],[0,WidtH]]
-    #m.add_region_from_polygon(outer_polygon, tags={'wall':[0,1,2], 'wave':[3]})
-
-    #domain.set_region(Set_region('mound', 'elevation', 100, location='unique vertices'))
 
     #Boundary
@@ -38 +31 @@
                             'wall',                           
                             'wall',
-                            'back',
+                            'back', # defines direchet boundary for back of test area
                             'wall',
                             'wall',
-                            'wave',
+                            'wave', # defines the wave boundary
                             '',                            
                             '']
@@ -51 +44 @@
     m.generate_mesh(maximum_triangle_area=maximum_triangle_area)
     triangle_count = m.get_triangle_count()
-  
-
     
     if mesh_file is None:    

development/momentum_sink/loop_buildings.py

@@ -10 +10 @@
 File_boundary, Dirichlet_boundary, Time_boundary, Transmissive_boundary
 from pyvolution.least_squares import fit_to_mesh_file, DEFAULT_ALPHA 
-
+import project
+from create_buildings import create_mesh
+from pmesh.mesh import importMeshFromFile
 import Numeric
 
@@ -19 +21 @@
 
 
-DR = [3,5,7,9]
+DR = [3,5,7,9,11,13,15,17]
 for depth in DR:
-    import project
     
-    from create_buildings import create_mesh
-    from pmesh.mesh import importMeshFromFile
     meshname = project.meshname
     outputname = project.outputname
     t0 = time.time()
-    meshname, triagle_count = cache(create_mesh,(100),
+    meshname, triagle_count = cache(create_mesh,(1000,depth),
                                     {'mesh_file':meshname,
                                      'triangles_in_name':True}
@@ -44 +43 @@
 
 
-    domain.set_name(project.basename + '_brd=%dm_%d' %(depth, triagle_count))
+    domain.set_name(project.basename + '_Orth(45)_D=%d_%d' %(depth, triagle_count))
     domain.set_datadir(project.outputdir)
     domain.store = True
-    domain.quantities_to_be_stored = ['stage']
+    domain.quantities_to_be_stored = ['stage', 'xmomentum', 'ymomentum']
 
     print 'Number of triangles = ', len(domain)
     print 'The extent is ', domain.get_extent()
+    print 'current building depth = ', depth
 
     #Setup Initial Conditions
@@ -66 +66 @@
     Bdb = Dirichlet_boundary([0,0,0])
     Bw = Time_boundary(domain=domain,
-                       f=lambda t: [(60<t<300)*4, 0, 0])
+                       f=lambda t: [(60<t<660)*4, 0, 0])
 
     domain.set_boundary( {'wall': Br,'wave': Bdw, 'back': Bdb, 'exterior':Bdw} ) 
@@ -73 +73 @@
     t0 = time.time()
 
-    for t in domain.evolve(yieldstep = 1, finaltime = 10):
+    for t in domain.evolve(yieldstep = 10, finaltime = 1000):
         domain.write_time()      
 

development/momentum_sink/loop_friction.py

@@ -24 +24 @@
 from pmesh.mesh import importMeshFromFile
 
-Friction = Numeric.arrayrange(50,80,1)
+#Friction = Numeric.arrayrange(0.025,1,0.025)
+Friction = [150,200,500,1000,2000,5000,10000]
 for fric in Friction:
     meshname = project_friction.meshname
@@ -56 +57 @@
                    )               
 
-
-    domain.set_name(project_friction.basename + '%d_%d' %(fric, triagle_count))
+    #fr=str(fric).strip('.')
+    domain.set_name(project_friction.basename + '%s_%d' %(str(fric).strip('.'), triagle_count))
     domain.set_datadir(project_friction.outputdir)
     domain.store = True

development/momentum_sink/run_buildings.py

@@ -15 +15 @@
 
 
-from create_buildings import create_mesh, depth
+from create_buildings import create_mesh
 #from building_generator import create_mesh
 
 from pmesh.mesh import importMeshFromFile
+
+depth = 17 # depth of building side to oncoming wave
+breadth = 17 # breadth of building, width of building to oncoming wave
+#from loop_buildings import depth
 
 meshname = project.meshname
@@ -25 +29 @@
 t0 = time.time()
 
-meshname, triagle_count = cache(create_mesh,(100),
+meshname, triagle_count = cache(create_mesh,(1000,depth),
                                 {'mesh_file':meshname,
                                  'triangles_in_name':True}
@@ -31 +35 @@
                                 #,evaluate = True      
                                 )
-#meshname = 'build.tsh'
+#meshname = 'test.tsh'
 #outputname = outputname[:-4] + '_' + str(triagle_count) + outputname[-4:]
 
@@ -43 +47 @@
                dependencies = [meshname]                    
                ,verbose = False
+               ,evaluate = True      
+
                )