Changeset 2429

Timestamp:

Feb 21, 2006, 8:33:21 AM (18 years ago)

Author:

nicholas

Message:

Updated CCS.py with successful implementation of 2nd normal test.

Location:

development/momentum_sink

Files:

• CCS.py (modified) (5 diffs)
• PCS.py (modified) (2 diffs)
• create_buildings.py (modified) (2 diffs)
• loop_friction.py (modified) (3 diffs)

development/momentum_sink/CCS.py

@@ -10 +10 @@
 # from run_friction import fric
 from pylab import *
-    
+from numerical_tools import norm, corr, err  
 
-swwfile_B = project.outputdir + sep  + 'Buildings_3662.sww'
+swwfile_B = project.outputdir + sep  + 'Buildings_3260.sww'
 
 
-gauge_depth = Numeric.arrayrange(0, 2.5*WidtH, 5)
+gauge_depth = Numeric.arrayrange(0, 1.5*WidtH, 20)
+print len(gauge_depth),"gauge_depth length"
 gauge_breadth = 100
 gauge_locations = []
@@ -27 +28 @@
 quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum']
 
-
-
 f_B = file_function(swwfile_B,
                   quantities = quantities,
@@ -34 +33 @@
                   verbose = True,
                   use_cache = True)
+
+Corr_OK=[]
 fric_OK=[]
-Friction = [1,2,3,5,7,10,20,30,50]
+least_norm=[1]
+norm_OK=[]
+Friction = Numeric.arrayrange(1,51,1)
 for fric in Friction:
     swwfile_F = project.outputdir + sep  + 'friction_n=' + repr(fric) + '_3135.sww'
@@ -48 +51 @@
     stages_B = []
     stages_F = []
-    elevations = []
-    momenta = []
-    velocity = []
     for i, GL in enumerate(gauge_depth):
 
@@ -70 +70 @@
         
         
-    
-    diff = abs(Numeric.array(stages_B) - Numeric.array(stages_F))
+    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)
+     
+    print " "
+    print norm_F, "norm N_F"
+    print " "
+    G = corr(Ar_F, Ar_B)
+    norm_OK.append(norm_F)
 
-    if max(diff) < 1.2:
-        print "_-_-_-_-_-_-_-_-_-_-_-_"
-        print fric, " this file is OK"
-        print "-_-_-_-_-_-_-_-_-_-_-_-"
-        fric_OK.append(fric)
-        #print swwfile_F, "this file OK"
+    if norm_F < least_norm:
+        least_norm = norm_F
+        fric_OK = fric
     else:
-        print fric, " this file is definelty NOT OK"
-print fric_OK, "THESE FRICTION VALUES ARE OK"
+        print fric, " NOT OK" 
 
     
+    
+print " "
+print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
+print fric_OK, " <- This Manning's n is OK"
+print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
+print least_norm, "2nd normal of Friction vs Buildings  "
+print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
 
-

development/momentum_sink/PCS.py

@@ -11 +11 @@
 from pylab import *
 
-fric = 1    
+fric = 40    
 #swwfile = project.newoutputname + '.sww'
 #swwfile = project.outputname
 
 #swwfile_B = project.outputdir + sep  + 'Buildings_3790.sww'
-#swwfile_B = project.outputdir + sep  + 'Buildings_3662.sww'
-swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
-swwfile_F = project.outputdir + sep  + 'friction_n=20_3135.sww'
+swwfile_B = project.outputdir + sep  + 'Buildings_3662.sww'
+#swwfile_B = project.outputdir + sep  + 'friction_n=10_3135.sww'
+swwfile_F = project.outputdir + sep  + 'friction_n=40_3135.sww'
 
-gauge_depth = Numeric.arrayrange(0, 2.5*WidtH, 5)
+gauge_depth = Numeric.arrayrange(0, 1.5*WidtH, 5)
 gauge_breadth = 100
 gauge_locations = []
@@ -45 +45 @@
 
 
-T = Numeric.arrayrange(0,500,20)
+T = Numeric.arrayrange(0,1000,20)
 #T = [ 20, 50, 100, 150, 200 ]
 

development/momentum_sink/create_buildings.py

@@ -8 +8 @@
 
 WidtH = 200 # width of boudary in metres
-depth = 20 # depth of building side to oncoming wave
-breadth = 20 # breadth of building
+depth = 15 # depth of building side to oncoming wave
+breadth = 15 # breadth of building
 print "building footprint"
 print depth * breadth , "m^2"
@@ -37 +37 @@
     # inner polygons => building boundaries
     
-    wh = depth/2 
-    
-    lh = breadth/2 
-          
+    whs = depth/2 
+    lhs = breadth/2 
+    #Th = (45 *(3.142/180)) # sets an initial rotation      
+    Th = 0
     ForDep = (0.2*WidtH) + (BL/2)
     RearDep = 1.2*WidtH 
         
-    Breadths = Numeric.arrayrange( (BL/2), WidtH, (BL))
+    Breadths = Numeric.arrayrange( -(BL/2), WidtH+(BL/2), (BL))
     print Breadths, "Breadths"
     Depths = Numeric.arrayrange( ForDep, RearDep, BL )
     print Depths, "Depths"
 
-    for D in Depths:
-        #Breadths = Breadths + BL/2 #Used to offset buildings
+    for i,D in enumerate(Depths):
+        Breadths = Breadths + ((-1)**i)*(BL/2) #Used to offset buildings
         for B in Breadths:
-            polygon = [[D-wh,B-lh],[D+wh,B-lh],[D+wh,B+lh],[D-wh,B+lh]]
+            wh1 = (-whs) * math.cos(Th) + (-lhs) * math.sin(Th)
+            lh1 = (-lhs) * math.cos(Th) - (-whs) * math.sin(Th)
+            wh2 = (+whs) * math.cos(Th) + (-lhs) * math.sin(Th)
+            lh2 = (-lhs) * math.cos(Th) - (+whs) * math.sin(Th)
+            wh3 = (+whs) * math.cos(Th) + (+lhs) * math.sin(Th)
+            lh3 = (+lhs) * math.cos(Th) - (+whs) * math.sin(Th)
+            wh4 = (-whs) * math.cos(Th) + (+lhs) * math.sin(Th)
+            lh4 = (+lhs) * math.cos(Th) - (-whs) * math.sin(Th)
+            polygon = [[D+wh1,B+lh1],[D+wh2,B+lh2],[D+wh3,B+lh3],[D+wh4,B+lh4]]
+            #polygon = [[D-wh,B-lh],[D+wh,B-lh],[D+wh,B+lh],[D-wh,B+lh]]
             m.add_hole_from_polygon(polygon, tags={'wall':[0,1,2,3]})# Adds holes with reflective boundaries.        
-        
+            Th = Th + (37.3 *(3.142/180)) # keeps rotating individual buildings.
 
             #print polygon

development/momentum_sink/loop_friction.py

@@ -7 +7 @@
 # set friction value for zone
 
-# fwojhrt
 import time
 
 import project_friction
-
+import Numeric
 from pyvolution.pmesh2domain import pmesh_to_domain_instance 
 from caching import cache
@@ -25 +24 @@
 from pmesh.mesh import importMeshFromFile
 
-Friction = [1,2,3,5,7,10,20,30,50]
+Friction = Numeric.arrayrange(1,50,1)
 for fric in Friction:
     meshname = project_friction.meshname
     outputname = project_friction.outputname
-
+    print " "
+    print "_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_"
+    print fric, "Manning friction value"
+    print "-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-"
     t0 = time.time()
 
@@ -84 +86 @@
     t0 = time.time()
 
-    for t in domain.evolve(yieldstep = 10, finaltime = 500):
+    for t in domain.evolve(yieldstep = 10, finaltime = 1000):
         domain.write_time()