Changeset 3190 for production/sydney_2006/run_sydney_smf.py
- Timestamp:
- Jun 21, 2006, 9:31:57 AM (18 years ago)
- File:
-
- 1 edited
-
production/sydney_2006/run_sydney_smf.py (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
production/sydney_2006/run_sydney_smf.py
r2640 r3190 27 27 from pyvolution.quantity import Quantity 28 28 from Numeric import allclose 29 from utilities.polygon import inside_polygon 29 30 30 31 # Application specific imports … … 70 71 71 72 # original issue to Benfield 72 #interior_res = 500073 #interior_regions = [[project.harbour_polygon_2, interior_res],74 #[project.botanybay_polygon_2, interior_res]]73 interior_res = 5000 74 interior_regions = [[project.harbour_polygon_2, interior_res], 75 [project.botanybay_polygon_2, interior_res]] 75 76 76 77 # used for finer mesh 77 interior_res1 = 5000 78 interior_res2 = 315 79 interior_regions = [[project.newpoly1, interior_res1], 80 [project.south1, interior_res1], 81 [project.finepolymanly, interior_res2], 82 [project.finepolyquay, interior_res2]] 83 84 print 'number of interior regions', len(interior_regions) 78 #interior_res1 = 5000 79 #interior_res2 = 315 80 #interior_regions = [[project.newpoly1, interior_res1], 81 # [project.south1, interior_res1], 82 # [project.finepolymanly, interior_res2], 83 # [project.finepolyquay, interior_res2]] 84 85 # used for coastal polygons constructed by GIS guys 86 def get_polygon_from_file(filename): 87 """ Function to read in output from GIS determined polygon 88 """ 89 fid = open(filename) 90 lines = fid.readlines() 91 fid.close() 92 93 polygon = [] 94 for line in lines[1:]: 95 fields = line.split(',') 96 x = float(fields[1]) 97 y = float(fields[2]) 98 polygon.append([x, y]) 99 100 return polygon 101 102 num_polygons = 9 103 fileext = '.csv' 104 filename = project.polygonptsfile 105 106 #interior_res = 1000 107 #interior_regions = [] 108 #bounding_polygon = project.diffpolygonall#project.demopoly 109 #count = 0 110 #for p in range(1, num_polygons+1): 111 # thefilename = filename + str(p) + fileext 112 # print 'reading in polygon points', thefilename 113 # interior_polygon = get_polygon_from_file(thefilename) 114 # interior_regions.append([interior_polygon, interior_res]) 115 # n = len(interior_polygon) 116 # # check interior polygon falls in bounding polygon 117 # if len(inside_polygon(interior_polygon, bounding_polygon, 118 # closed = True, verbose = False)) <> len(interior_polygon): 119 # print 'WARNING: interior polygon %d is outside bounding polygon' %(p) 120 # count += 1 121 # # check for duplicate points in interior polygon 122 123 print 'number of interior polygons: ', len(interior_regions) 124 #if count == 0: print 'interior polygons OK' 85 125 86 126 #FIXME: Fix caching of this one once the mesh_interface is ready … … 89 129 # original + refined region 90 130 _ = cache(create_mesh_from_regions, 91 project.diffpolygonall, 131 #project.demopoly, 132 project.diffpolygonall2, 133 #{'boundary_tags': {'bottom': [0], 134 # 'right1': [1], 'right0': [2], 135 # 'right2': [3], 'top': [4], 'left1': [5], 136 # 'left2': [6], 'left3': [7]}, 92 137 {'boundary_tags': {'bottom': [0], 93 ' right1': [1], 'right0': [2],94 ' right2': [3], 'top': [4], 'left1': [5],138 'bottom1': [1], 'right': [2], 139 'top1': [3], 'top': [4], 'left1': [5], 95 140 'left2': [6], 'left3': [7]}, 141 #{'boundary_tags': {'bottom': [0], 'right': [1], 142 # 'top': [2], 'left': [3]}, 96 143 'maximum_triangle_area': 100000, 97 144 'filename': meshname, … … 126 173 radius=3330, 127 174 dphi=0.23, 128 x0=project.slump_origin [0],129 y0=project.slump_origin [1],175 x0=project.slump_origin2[0], 176 y0=project.slump_origin2[1], 130 177 alpha=0.0, 131 domain=domain )132 178 domain=domain, 179 verbose=True) 133 180 134 181 #------------------------------------------------------------------------------- … … 138 185 # apply slump after 30 mins, initialise to water level of tide = 0 139 186 domain.set_quantity('stage', 0.0) 140 domain.set_quantity('friction', 0.0 3)187 domain.set_quantity('friction', 0.04) 141 188 domain.set_quantity('elevation', 142 189 filename = project.combineddemname + '.pts', … … 158 205 # 'right2': Br, 'top': Br, 'left1': Br, 159 206 # 'left2': Br, 'left3': Br} ) 160 207 # for new tests 4 April 2006 208 #domain.set_boundary( {'bottom': Br, 'bottom1': Br, 'right': Br, 209 # 'top1': Br, 'top': Br, 'left1': Br, 210 # 'left2': Br, 'left3': Br} ) 211 161 212 # enforce Dirichlet BC - from 30/03/06 Benfield visit 162 domain.set_boundary( {'bottom': Bd, ' right1': Bd, 'right0': Bd,163 ' right2': Bd, 'top': Bd, 'left1': Bd,213 domain.set_boundary( {'bottom': Bd, 'bottom1': Bd, 'right': Bd, 214 'top1': Bd, 'top': Bd, 'left1': Bd, 164 215 'left2': Bd, 'left3': Bd} ) 165 216 166 217 # increasingly finer interior regions 167 #domain.set_boundary( {'bottom': B r, 'right': Br, 'left': Br, 'top': Br} )218 #domain.set_boundary( {'bottom': Bd, 'right': Bd, 'left': Bd, 'top': Bd} ) 168 219 169 220 … … 177 228 fid = open(thisfile, 'w') 178 229 179 for t in domain.evolve(yieldstep = 120, finaltime = 18000): 230 # save every 10 secs leading up to slump initiation 231 for t in domain.evolve(yieldstep = 10, finaltime = 60): # 6 steps 180 232 domain.write_time() 181 233 domain.write_boundary_statistics(tags = 'bottom') 182 183 # calculate integral 184 thisstagestep = domain.get_quantity('stage') 185 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 186 fid.write(s) 187 188 # add slump after 30 mins 189 if allclose(t, 30*60): 234 # calculate integral 235 thisstagestep = domain.get_quantity('stage') 236 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 237 fid.write(s) 238 # add slump after 30 secs 239 if allclose(t, 30): 190 240 slump = Quantity(domain) 191 241 slump.set_values(tsunami_source) 192 242 domain.set_quantity('stage', slump + thisstagestep) 193 243 #test_stage = domain.get_quantity('stage') 244 #test_elevation = domain.get_quantity('elevation') 245 #test_depth = test_stage - test_elevation 246 #test_max = max(test_depth.get_values()) 247 #print 'testing', test_max 248 249 import sys; sys.exit() 250 251 # save every two minutes leading up to interesting period 252 for t in domain.evolve(yieldstep = 120, finaltime = 660, # steps 253 skip_initial_step = True): 254 domain.write_time() 255 domain.write_boundary_statistics(tags = 'bottom') 256 # calculate integral 257 thisstagestep = domain.get_quantity('stage') 258 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 259 fid.write(s) 260 261 262 # save every thirty secs during interesting period 263 for t in domain.evolve(yieldstep = 60, finaltime = 5000, # steps 264 skip_initial_step = True): 265 domain.write_time() 266 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage') 267 # calculate integral 268 thisstagestep = domain.get_quantity('stage') 269 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 270 fid.write(s) 271 272 273 import sys; sys.exit() 274 # save every two mins for next 5000 secs 275 for t in domain.evolve(yieldstep = 120, finaltime = 10000, # about 42 steps 276 skip_initial_step = True): 277 domain.write_time() 278 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage') 279 # calculate integral 280 thisstagestep = domain.get_quantity('stage') 281 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 282 fid.write(s) 283 284 # save every half hour to end of simulation 285 for t in domain.evolve(yieldstep = 1800, finaltime = 10*60*60, # 14 steps 286 skip_initial_step = True): 287 domain.write_time() 288 domain.write_boundary_statistics(tags = 'bottom') #quantities = 'stage' 289 # calculate integral 290 thisstagestep = domain.get_quantity('stage') 291 s = '%.2f, %.2f\n' %(t, thisstagestep.get_integral()) 292 fid.write(s) 194 293 195 294 print 'That took %.2f seconds' %(time.time()-t0)
Note: See TracChangeset
for help on using the changeset viewer.
