Changeset 5038

Timestamp:

Feb 18, 2008, 1:32:33 PM (17 years ago)

Author:

nick

Message:

updates of run_circular.py

File:

• anuga_validation/circular_island/run_circular.py (modified) (14 diffs)

anuga_validation/circular_island/run_circular.py

@@ -2 +2 @@
 
 This script sets up 2D circular island.....
+
+use 'res' parameter of 1 to make run fast.
 
 """
@@ -14 +16 @@
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
 from anuga.utilities.polygon import read_polygon#, plot_polygons
-#from cmath import cos,sin,cosh,pi
 from math import cos,pi,sin,tan,sqrt
 from Numeric import array, zeros, Float, allclose,resize
@@ -94 +95 @@
           247.5, 270.0, 292.5, 315.0, 337.5)
 
-#x = 12.96
-#y = 13.80
 center_x = 12.96
 center_y = 13.80
 
-#wavelenght = 
-
-r1 = 1.5
-r2 = 2.5
 d=0.75
-#res=.05
+res=0.05
 res=1.
-L=4.11
-
+#L=3.6 and therefore the gauges 1,2,3 and 4 are located 1.8 from the toe of the island
+# 13.8 - 3.6(radii of island) - 1.8 = 8.4
+L=8.4
 
 #create polygons (circles) to have higher resolution
@@ -117 +113 @@
     angle = ((angle-180)/180.0)*pi
 
-    p = get_xy(center_x,center_y,3.6,angle)
+    p = get_xy(center_x,center_y,4.0,angle)
     poly1.append([p[0],p[1]])
 
@@ -127 +123 @@
     #convert Degs to Rads
     angle = ((angle-180)/180.0)*pi
-    p = get_xy(center_x,center_y,2.4,angle)
+    p = get_xy(center_x,center_y,2.6,angle)
     poly2.append([p[0],p[1]])
     
@@ -137 +133 @@
     #convert Degs to Rads
     angle = ((angle-180)/180.0)*pi
-    p = get_xy(center_x,center_y,1.6,angle)
+#    p = get_xy(center_x,center_y,1.6,angle)
+    p = get_xy(center_x,center_y,1.5,angle)
     poly3.append([p[0],p[1]])
     
@@ -143 +140 @@
 internal_poly.append([poly3,1*res])
 
-poly4 = []
-for i,angle in enumerate(angles1):
-    #convert Degs to Rads
-    angle = ((angle-180)/180.0)*pi
-    p = get_xy(center_x,center_y,1.1,angle)
-    poly4.append([p[0],p[1]])
-    
-poly.append(poly4)
-internal_poly.append([poly2,.01*res])
-
 #plot_polygons(poly,'poly.png')
 
-#poly_all= [[0,7],[0,25],[30,25],[30,7]]
 poly_all= [[0,L],[0,25],[30,25],[30,L]]
 create_mesh_from_regions(poly_all,
                          boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]},
-                         maximum_triangle_area=1,
+                         maximum_triangle_area=1*res,
                          interior_regions=internal_poly,
                          filename='temp.msh',
@@ -180 +166 @@
 #print attribute_dic['number'][10]
 
-def test(x,y):
-#    center_x = 12.96
-#    center_y = 13.80
-#    center_y = 6.80
-
-    list_xy = (center_x-x)**2+(center_y-L-y)**2
+def circular_island_elevation(x,y):
+
+    list_xy_square = (center_x-x)**2+(center_y-L-y)**2
     print 'x',min(x),max(x)
     print 'y',min(y),max(y)
-    z1=[]
-    for xy in list_xy:
+    list_z=[]
+    for xy in list_xy_square:
         # this finds the int that relates to the correct sqrt in the table
         nn=int(round(xy,2)*100)
-        zz=-float(attribute_dic['sqrt'][nn])
+        #the square root is to find the distance from the center to the current xy
+        distance = float(attribute_dic['sqrt'][nn])
+#        zz=-float(attribute_dic['sqrt'][nn])
 #        print nn,zz
-
-#    z = -square_root((center_x-x)**2+(center_y-y)**2)
-        z = zz*.25+0.8975
+        
+        # slope of cone is 1/4 and has a radii of 3.6 therefore the cone is 0.9 high. 
+        # So to have the elevation positive 0.90 has been added
+#        z = zz*.25+0.90
+        z = -distance*.25+0.90
+#        z = zz*.25+0.8975
 
         if z <0:
@@ -203 +191 @@
             z=0.625
 #        print 'hello',z
-        z1.append(z)
-    return z1
-
-domain.set_quantity('elevation',test, verbose=True)
+        list_z.append(z)
+    return list_z
+
+domain.set_quantity('elevation', circular_island_elevation, verbose=True)
 
 
@@ -227 +215 @@
 
 # Create boundary function from timeseries provided in file
-#function = file_function(project.boundary_filename,
-#                         domain, verbose=True)
-
-# Create and assign boundary objects
+
 boundary_filename="ts2cnew1_input"
 prepare_timeboundary(boundary_filename+'.txt')
@@ -245 +230 @@
 Bd = Dirichlet_boundary([water_height,0,0])
 domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd})
-domain.starttime = 10
+#domain.set_time(20.0)
 #-------------------------
 # Evolve through time
@@ -273 +258 @@
           95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0,
           247.5, 270.0, 292.5, 315.0, 337.5)
-#x = 12.96
-#y = 13.80
-r1 = 1.1
-r2 = 3.6
+r1 = 1.5
+r2 = 2.6
+#degrees to check
 d=2
 
@@ -283 +267 @@
 for i,angle in enumerate(angles):
     angle = ((angle-180)/180.0)*pi
-#    angle = angle1*3.14157
     p1 = get_xy(center_x,center_y,r1,angle-0.01745*d)
     p2 = get_xy(center_x,center_y,r2,angle-0.01745*d)
@@ -295 +278 @@
     run_up, location = get_maximum_inundation_data(filename=model_name+'.sww',polygon=poly_segment, verbose=False)
     
-    print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, run_up-water_height*100, location[0],location[1]) 
-
-
-
-
-
-
-
+    print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, (run_up-water_height)*100, location[0],location[1]) 
+
+
+
+
+
+
+