Changeset 5977
- Timestamp:
- Nov 19, 2008, 3:24:59 PM (16 years ago)
- Location:
- anuga_work/development/ted_rigby
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/ted_rigby/model_data.py
r5973 r5977 322 322 # 323 323 # OUTPUTS OF THIS CODE BLOCK ARE; 324 # mat_roughness_data _list[[RID],[n0] etc.] -- listof properties for each surface (RID)325 # mat_RID_list[RID] 326 # mat_poly_list[[x0,y0],[x1,y1]...] 324 # mat_roughness_data[RID] -- dictionary of properties for each surface (RID) 325 # mat_RID_list[RID] -- surface roughness Ids (RID) for each poly 326 # mat_poly_list[[x0,y0],[x1,y1]...] -- list of vertices for each roughness poly 327 327 #-------------------------------------------------------------------------------- 328 328 … … 330 330 print basename+ ' >>>> Assemble the surface roughness data from Tuflows %s and %s ' %(tmf_filename, mat_filename) 331 331 332 # Get the tuflow roughness data for each surface class 332 333 mat_roughness_data_list=[] 333 mat_roughness_data_list=get_tuflow_data.get_tuflow_tmf(tmf_filename,verbose=model_verbose) # get the tuflow roughness data for each surface class 334 334 mat_roughness_data_list=get_tuflow_data.get_tuflow_tmf(tmf_filename,verbose=model_verbose) 335 336 # Turn material roughness data list into a dictionary 337 mat_roughness_data = {} 338 for entry in mat_roughness_data_list: 339 key = entry[0] # rid 340 mat_roughness_data[key] = [entry[1], # n0 341 entry[4], # d1 342 entry[5], # n1 343 entry[6], # d2 344 entry[7]] # n2 345 346 347 348 335 349 mat_RID_list=[] # Int list of rough poly IDs-one for each poly 336 350 mat_poly_list=[] # list of float list of poly coords -
anuga_work/development/ted_rigby/model_run.py
r5975 r5977 24 24 import time 25 25 import sys 26 from Numeric import zeros, Float 26 27 27 28 # Related ANUGA modules … … 35 36 from anuga.abstract_2d_finite_volumes.quantity import Quantity 36 37 from anuga.pmesh.mesh_interface import create_mesh_from_regions 37 from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function 38 from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function, inside_polygon 38 39 from anuga.alpha_shape.alpha_shape import alpha_shape_via_files 39 40 … … 76 77 77 78 78 domain79 domain.forcing_terms.append(hydrograph)80 81 print model_data.basename+' >>>> Completed assignment of %i inflow hydrographs at t= %.2f hours ' %(len(model_data.infl82 79 print model_data.basename+' >>>> Completed mesh creation at t= %.2f hours ' %(float(time.time()-t0)/3600.0) 83 80 … … 245 242 #-------------------------------------------------------------------------------------- 246 243 247 print model_data.basename+' >>>> Creating and assigning mat_RID (surface Roughness ID)instance for use in evolve' 248 249 # Create an instance of a new domain quantity called mat_RID to store the surface Roughness IDs - fill initially with zeros 250 mat_RID = Quantity(domain) 251 if model_data.model_verbose : 252 print ' ++++ Created a zeroed domain quantity of %s and length %i called mat_RID' %( type(mat_RID),len(mat_RID) ) 253 254 # Set all of the RID values initially to one fixed default value (makes sure all elements have a RID) 255 mat_RID.set_values(numeric=model_data.default_RID,location='centroids',verbose=True) 256 if model_data.model_verbose : 257 print ' ++++ Initially filled mat_RID with the specified default RID value of %i ' %( model_data.default_RID ) 258 259 # Assign model_data.mat_RID_list values to domain.from model_data.mat_poly_list 260 for i in range(len(model_data.mat_RID_list)) : 261 mat_RID.set_values(numeric=model_data.mat_RID_list[i],polygon=model_data.mat_poly_list[i]) 262 if model_data.model_verbose : 263 print ' ++++ Upgraded mat_RID spatially with the values set in the mif/mid files ' 264 265 266 # Dump some RIDS to confirm assignment occured as intended 267 if model_data.model_verbose : 268 print ' ++++ Check printout of RID assignments to centroids of elements nos 100 200 and 300' 269 print ' ++++ Element 100s RIDs ', mat_RID.get_values(location='centroids',indices=[100]) 270 print ' ++++ Element 200s RIDs ', mat_RID.get_values(location='centroids',indices=[200]) 271 print ' ++++ Element 300s RIDs ', mat_RID.get_values(location='centroids',indices=[300]) 272 273 print model_data.basename+' >>>> Completed assignment of material (surface roughness) RIDs and default n to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0) 244 print model_data.basename+' >>>> Creating and assigning surface roughness data for use in evolve' 245 246 N = len(domain) 247 material_variables = zeros((N, 5), Float) 248 249 # Set default material variables n0, d1, n1, d2, n2 250 material_variables[:,0] = 0.035 # n0 251 material_variables[:,1] = 0.3 # d1 252 material_variables[:,2] = 0.050 # n1 253 material_variables[:,3] = 1.5 # d2 254 material_variables[:,4] = 0.030 # n2 255 256 257 # Assign polygonal data to material_variables from model_data.mat_poly_list 258 259 260 points = domain.get_centroid_coordinates(absolute=True) 261 262 for k, poly in enumerate(model_data.mat_poly_list): 263 264 # Get indices of triangles inside polygon k 265 indices = inside_polygon(points, poly, closed=True) 266 267 # Index into material data for polygon k 268 rid = model_data.mat_RID_list[k] 269 270 # Material data for polygon 271 n0, d1, n1, d2, n2 = model_data.mat_roughness_data[rid] 272 273 for i in indices: 274 material_variables[i, 0] = n0 275 material_variables[i, 1] = d1 276 material_variables[i, 2] = n1 277 material_variables[i, 3] = d2 278 material_variables[i, 4] = n2 279 280 281 282 print ' ++++ Upgraded material_variables spatially with the material values set in the mif/mid files ' 283 284 285 print model_data.basename+' >>>> Completed assignment of material (surface roughness) and default n to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0) 274 286 275 287 print model_data.basename+' >>>> Completed domain construction at t= %.2f hours ' %(float(time.time()-t0)/3600.0) 288 276 289 277 290 … … 316 329 print ' ++++ Recomputing friction at t= ', t 317 330 318 noelementsupdated = 0 # intialise update check 331 319 332 # rebuild the 'friction' values adjusted for depth 320 333 for i in domain.get_wet_elements(): # loop for each wet element in domain 321 334 322 element_RID_value=mat_RID.get_values(location='centroids', 323 indices=[i]) # get current elements centroidal RID value 324 325 # print ' ++++ Element %i centroidal element_RID_value is %s' %( i,element_RID_value[0] ) 326 for k in range(len(model_data.mat_roughness_data_list)): # loop for each matl in tmf list 327 if int(element_RID_value[0]) == int( model_data.mat_roughness_data_list[k][0] ) : # got tmf line for that matl 328 # print ' ++++ found tmf match for RID %i of %s ' %( int(element_RID_value[0]),model_data.mat_roughness_data_list[k][8].rstrip(' ') ) 329 d1 = float(model_data.mat_roughness_data_list[k][4]) 330 n1 = float(model_data.mat_roughness_data_list[k][5]) 331 d2 = float(model_data.mat_roughness_data_list[k][6]) 332 n2 = float(model_data.mat_roughness_data_list[k][7]) 333 noelementsupdated = noelementsupdated + 1 # increment update counter start at zeroth element 334 if depth.get_values(location='centroids', 335 indices=[i]) <= d1 : # recompute friction values from depth for this element get_tuflow_data.py 336 337 domain.set_quantity('friction', 338 numeric=n1, 339 location='centroids', 340 indices=[i], 341 verbose=model_data.anuga_verbose) 342 343 elif depth.get_values(location='centroids', 344 indices=[i]) >= d2 : 345 domain.set_quantity('friction',numeric=n2,location='centroids', 346 indices=[i],verbose=model_data.anuga_verbose) 347 else : 348 n3 = n1+(n2-n1)/(d2-d1)* ( depth.get_values(location='centroids',indices=[i])[0] ) 349 domain.set_quantity('friction',numeric=n3,location='centroids', 350 indices=[i],verbose=model_data.anuga_verbose) 351 335 # Get roughness variables 336 n0 = material_variables[i,0] 337 d1 = material_variables[i,1] 338 n1 = material_variables[i,2] 339 d2 = material_variables[i,3] 340 n2 = material_variables[i,4] 341 342 343 # Recompute friction values from depth for this element get_tuflow_data.py 344 d = depth.get_values(location='centroids', indices=[i])[0] 345 346 if d <= d1: 347 depth_dependent_friction = n1 348 elif d >= d2: 349 depth_dependent_friction = n2 350 else: 351 depth_dependent_friction = n1+(n2-n1)/(d2-d1)*d 352 353 354 domain.set_quantity('friction', 355 numeric=depth_dependent_friction, 356 location='centroids', 357 indices=[i], 358 verbose=model_data.anuga_verbose) 359 352 360 353 361 if model_data.model_verbose : 354 355 362 friction = domain.get_quantity('friction') 356 363 357 print " ++++ Check some RID,depth, friction values returned during this timestep " 358 print " element 100 ", mat_RID.get_values(location='centroids',indices=[100]), depth.get_values(location='centroids',indices=[100]), friction.get_values(location='centroids',indices=[100]) 359 360 print " element 200 ", mat_RID.get_values(location='centroids',indices=[200]), depth.get_values(location='centroids',indices=[200]), friction.get_values(location='centroids',indices=[200]) 361 362 print " element 300 ", mat_RID.get_values(location='centroids',indices=[300]), depth.get_values(location='centroids',indices=[300]), friction.get_values(location='centroids',indices=[300] ) 363 364 365 364 # Print something 365 366 366 367 367
Note: See TracChangeset
for help on using the changeset viewer.