Changeset 5977


Ignore:
Timestamp:
Nov 19, 2008, 3:24:59 PM (14 years ago)
Author:
ole
Message:

Running version of variable friction

Location:
anuga_work/development/ted_rigby
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_work/development/ted_rigby/model_data.py

    r5973 r5977  
    322322#
    323323# OUTPUTS OF THIS CODE BLOCK ARE;
    324 # mat_roughness_data_list[[RID],[n0] etc.] -- list 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
     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
    327327#--------------------------------------------------------------------------------
    328328
     
    330330print basename+ ' >>>> Assemble the surface roughness data from Tuflows %s and %s ' %(tmf_filename, mat_filename)
    331331
     332# Get the tuflow roughness data for each surface class
    332333mat_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 
     334mat_roughness_data_list=get_tuflow_data.get_tuflow_tmf(tmf_filename,verbose=model_verbose)
     335
     336# Turn material roughness data list into a dictionary
     337mat_roughness_data = {}
     338for 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   
    335349mat_RID_list=[]                                                         # Int list of rough poly IDs-one for each poly
    336350mat_poly_list=[]                                                        # list of float list of poly coords
  • anuga_work/development/ted_rigby/model_run.py

    r5975 r5977  
    2424import time
    2525import sys
     26from Numeric import zeros, Float
    2627
    2728# Related ANUGA modules
     
    3536from anuga.abstract_2d_finite_volumes.quantity import Quantity
    3637from anuga.pmesh.mesh_interface import create_mesh_from_regions
    37 from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function
     38from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function, inside_polygon
    3839from anuga.alpha_shape.alpha_shape import alpha_shape_via_files
    3940
     
    7677
    7778
    78 domain
    79     domain.forcing_terms.append(hydrograph)
    80 
    81 print model_data.basename+' >>>> Completed assignment of %i inflow  hydrographs at t= %.2f hours ' %(len(model_data.infl
    8279print model_data.basename+' >>>> Completed mesh creation at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
    8380
     
    245242#--------------------------------------------------------------------------------------
    246243
    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)
     244print model_data.basename+' >>>> Creating and assigning surface roughness data for use in evolve'
     245
     246N = len(domain)
     247material_variables = zeros((N, 5), Float)
     248
     249# Set default material variables n0, d1, n1, d2, n2
     250material_variables[:,0] = 0.035 # n0
     251material_variables[:,1] = 0.3   # d1 
     252material_variables[:,2] = 0.050 # n1 
     253material_variables[:,3] = 1.5   # d2 
     254material_variables[:,4] = 0.030 # n2 
     255
     256
     257# Assign polygonal data to material_variables from model_data.mat_poly_list
     258
     259
     260points = domain.get_centroid_coordinates(absolute=True)
     261
     262for 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
     282print '         ++++ Upgraded material_variables spatially with the material values set in the mif/mid files '
     283
     284   
     285print model_data.basename+' >>>> Completed assignment of material (surface roughness) and default n to domain at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
    274286
    275287print model_data.basename+' >>>> Completed domain construction at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
     288
    276289
    277290
     
    316329    print '         ++++ Recomputing friction at t= ', t
    317330       
    318     noelementsupdated = 0                                                      # intialise update check
     331
    319332    # rebuild the 'friction' values adjusted for depth
    320333    for i in domain.get_wet_elements():                                        # loop for each wet element in domain
    321334   
    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             
    352360             
    353361    if model_data.model_verbose :
    354    
    355362        friction = domain.get_quantity('friction')
    356363       
    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       
    366366       
    367367
Note: See TracChangeset for help on using the changeset viewer.