Changeset 5038
 Timestamp:
 Feb 18, 2008, 1:32:33 PM (17 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_validation/circular_island/run_circular.py
r5036 r5038 2 2 3 3 This script sets up 2D circular island..... 4 5 use 'res' parameter of 1 to make run fast. 4 6 5 7 """ … … 14 16 from anuga.pmesh.mesh_interface import create_mesh_from_regions 15 17 from anuga.utilities.polygon import read_polygon#, plot_polygons 16 #from cmath import cos,sin,cosh,pi17 18 from math import cos,pi,sin,tan,sqrt 18 19 from Numeric import array, zeros, Float, allclose,resize … … 94 95 247.5, 270.0, 292.5, 315.0, 337.5) 95 96 96 #x = 12.9697 #y = 13.8098 97 center_x = 12.96 99 98 center_y = 13.80 100 99 101 #wavelenght =102 103 r1 = 1.5104 r2 = 2.5105 100 d=0.75 106 #res=.05101 res=0.05 107 102 res=1. 108 L=4.11 109 103 #L=3.6 and therefore the gauges 1,2,3 and 4 are located 1.8 from the toe of the island 104 # 13.8  3.6(radii of island)  1.8 = 8.4 105 L=8.4 110 106 111 107 #create polygons (circles) to have higher resolution … … 117 113 angle = ((angle180)/180.0)*pi 118 114 119 p = get_xy(center_x,center_y, 3.6,angle)115 p = get_xy(center_x,center_y,4.0,angle) 120 116 poly1.append([p[0],p[1]]) 121 117 … … 127 123 #convert Degs to Rads 128 124 angle = ((angle180)/180.0)*pi 129 p = get_xy(center_x,center_y,2. 4,angle)125 p = get_xy(center_x,center_y,2.6,angle) 130 126 poly2.append([p[0],p[1]]) 131 127 … … 137 133 #convert Degs to Rads 138 134 angle = ((angle180)/180.0)*pi 139 p = get_xy(center_x,center_y,1.6,angle) 135 # p = get_xy(center_x,center_y,1.6,angle) 136 p = get_xy(center_x,center_y,1.5,angle) 140 137 poly3.append([p[0],p[1]]) 141 138 … … 143 140 internal_poly.append([poly3,1*res]) 144 141 145 poly4 = []146 for i,angle in enumerate(angles1):147 #convert Degs to Rads148 angle = ((angle180)/180.0)*pi149 p = get_xy(center_x,center_y,1.1,angle)150 poly4.append([p[0],p[1]])151 152 poly.append(poly4)153 internal_poly.append([poly2,.01*res])154 155 142 #plot_polygons(poly,'poly.png') 156 143 157 #poly_all= [[0,7],[0,25],[30,25],[30,7]]158 144 poly_all= [[0,L],[0,25],[30,25],[30,L]] 159 145 create_mesh_from_regions(poly_all, 160 146 boundary_tags={'wall': [0,2],'wave': [3], 'buffer':[1]}, 161 maximum_triangle_area=1 ,147 maximum_triangle_area=1*res, 162 148 interior_regions=internal_poly, 163 149 filename='temp.msh', … … 180 166 #print attribute_dic['number'][10] 181 167 182 def test(x,y): 183 # center_x = 12.96 184 # center_y = 13.80 185 # center_y = 6.80 186 187 list_xy = (center_xx)**2+(center_yLy)**2 168 def circular_island_elevation(x,y): 169 170 list_xy_square = (center_xx)**2+(center_yLy)**2 188 171 print 'x',min(x),max(x) 189 172 print 'y',min(y),max(y) 190 z1=[]191 for xy in list_xy :173 list_z=[] 174 for xy in list_xy_square: 192 175 # this finds the int that relates to the correct sqrt in the table 193 176 nn=int(round(xy,2)*100) 194 zz=float(attribute_dic['sqrt'][nn]) 177 #the square root is to find the distance from the center to the current xy 178 distance = float(attribute_dic['sqrt'][nn]) 179 # zz=float(attribute_dic['sqrt'][nn]) 195 180 # print nn,zz 196 197 # z = square_root((center_xx)**2+(center_yy)**2) 198 z = zz*.25+0.8975 181 182 # slope of cone is 1/4 and has a radii of 3.6 therefore the cone is 0.9 high. 183 # So to have the elevation positive 0.90 has been added 184 # z = zz*.25+0.90 185 z = distance*.25+0.90 186 # z = zz*.25+0.8975 199 187 200 188 if z <0: … … 203 191 z=0.625 204 192 # print 'hello',z 205 z1.append(z)206 return z1207 208 domain.set_quantity('elevation', test, verbose=True)193 list_z.append(z) 194 return list_z 195 196 domain.set_quantity('elevation', circular_island_elevation, verbose=True) 209 197 210 198 … … 227 215 228 216 # Create boundary function from timeseries provided in file 229 #function = file_function(project.boundary_filename, 230 # domain, verbose=True) 231 232 # Create and assign boundary objects 217 233 218 boundary_filename="ts2cnew1_input" 234 219 prepare_timeboundary(boundary_filename+'.txt') … … 245 230 Bd = Dirichlet_boundary([water_height,0,0]) 246 231 domain.set_boundary({'wave': Bts, 'wall': Br, 'buffer':Bd}) 247 domain.starttime = 10 232 #domain.set_time(20.0) 248 233 # 249 234 # Evolve through time … … 273 258 95.0, 100.0, 105.0, 112.5, 135.0, 157.5, 180.0, 202.5, 225.0, 274 259 247.5, 270.0, 292.5, 315.0, 337.5) 275 #x = 12.96 276 #y = 13.80 277 r1 = 1.1 278 r2 = 3.6 260 r1 = 1.5 261 r2 = 2.6 262 #degrees to check 279 263 d=2 280 264 … … 283 267 for i,angle in enumerate(angles): 284 268 angle = ((angle180)/180.0)*pi 285 # angle = angle1*3.14157286 269 p1 = get_xy(center_x,center_y,r1,angle0.01745*d) 287 270 p2 = get_xy(center_x,center_y,r2,angle0.01745*d) … … 295 278 run_up, location = get_maximum_inundation_data(filename=model_name+'.sww',polygon=poly_segment, verbose=False) 296 279 297 print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, run_upwater_height*100, location[0],location[1])298 299 300 301 302 303 304 280 print 'maxd %.3f, %.6f, %.6f, %.5s, %.5s '%(((angle/pi)*180)+180, run_up*100, (run_upwater_height)*100, location[0],location[1]) 281 282 283 284 285 286 287
Note: See TracChangeset
for help on using the changeset viewer.