Changeset 5977

Timestamp:

Nov 19, 2008, 3:24:59 PM (15 years ago)

Author:

ole

Message:

Running version of variable friction

Location:

anuga_work/development/ted_rigby

Files:

• model_data.py (modified) (2 diffs)
• model_run.py (modified) (5 diffs)

anuga_work/development/ted_rigby/model_data.py

@@ -322 +322 @@
 #
 # OUTPUTS OF THIS CODE BLOCK ARE;
-# mat_roughness_data_list[[RID],[n0] etc.] -- list of properties for each surface (RID)
-# mat_RID_list[RID]                        --  surface roughness Ids (RID) for each poly 
-# mat_poly_list[[x0,y0],[x1,y1]...]        --  list of vertices for each roughness poly
+# mat_roughness_data[RID]            --  dictionary of properties for each surface (RID)
+# mat_RID_list[RID]                  --  surface roughness Ids (RID) for each poly 
+# mat_poly_list[[x0,y0],[x1,y1]...]  --  list of vertices for each roughness poly
 #--------------------------------------------------------------------------------
 
@@ -330 +330 @@
 print basename+ ' >>>> Assemble the surface roughness data from Tuflows %s and %s ' %(tmf_filename, mat_filename)
 
+# Get the tuflow roughness data for each surface class
 mat_roughness_data_list=[]
-mat_roughness_data_list=get_tuflow_data.get_tuflow_tmf(tmf_filename,verbose=model_verbose)    # get the tuflow roughness data for each surface class
-
+mat_roughness_data_list=get_tuflow_data.get_tuflow_tmf(tmf_filename,verbose=model_verbose)
+
+# Turn material roughness data list into a dictionary
+mat_roughness_data = {}
+for entry in mat_roughness_data_list:
+    key = entry[0]                        # rid 
+    mat_roughness_data[key] = [entry[1],  # n0 
+                               entry[4],  # d1
+                               entry[5],  # n1
+                               entry[6],  # d2
+                               entry[7]]  # n2
+                               
+    
+
+    
 mat_RID_list=[]                                                         # Int list of rough poly IDs-one for each poly
 mat_poly_list=[]                                                        # list of float list of poly coords

anuga_work/development/ted_rigby/model_run.py

@@ -24 +24 @@
 import time
 import sys
+from Numeric import zeros, Float
 
 # Related ANUGA modules
@@ -35 +36 @@
 from anuga.abstract_2d_finite_volumes.quantity import Quantity
 from anuga.pmesh.mesh_interface import create_mesh_from_regions
-from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function
+from anuga.utilities.polygon import read_polygon, plot_polygons, Polygon_function, inside_polygon
 from anuga.alpha_shape.alpha_shape import alpha_shape_via_files
 
@@ -76 +77 @@
 
 
-domain
-    domain.forcing_terms.append(hydrograph)
-
-print model_data.basename+' >>>> Completed assignment of %i inflow  hydrographs at t= %.2f hours ' %(len(model_data.infl
 print model_data.basename+' >>>> Completed mesh creation at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
 
@@ -245 +242 @@
 #--------------------------------------------------------------------------------------
 
-print model_data.basename+' >>>> Creating and assigning mat_RID (surface Roughness ID)instance for use in evolve'
-
-# Create an instance of a new domain quantity called mat_RID to store the surface Roughness IDs - fill initially with zeros
-mat_RID = Quantity(domain)
-if model_data.model_verbose :
-    print '         ++++ Created a zeroed domain quantity of  %s and length %i called mat_RID' %( type(mat_RID),len(mat_RID) )
-
-# Set  all of the RID values initially to one fixed default value (makes sure all elements have a RID)
-mat_RID.set_values(numeric=model_data.default_RID,location='centroids',verbose=True)
-if model_data.model_verbose :
-    print '         ++++ Initially filled mat_RID with the specified default RID value of %i ' %( model_data.default_RID )
-
-# Assign  model_data.mat_RID_list values to domain.from model_data.mat_poly_list
-for i in range(len(model_data.mat_RID_list)) :
-    mat_RID.set_values(numeric=model_data.mat_RID_list[i],polygon=model_data.mat_poly_list[i])
-if model_data.model_verbose :
-    print '         ++++ Upgraded mat_RID spatially with the values set in the mif/mid files ' 
-
-
-# Dump some RIDS to confirm assignment occured as intended
-if model_data.model_verbose :
-    print '         ++++ Check printout of RID assignments to centroids of elements nos 100 200 and 300'
-    print '         ++++ Element 100s RIDs ', mat_RID.get_values(location='centroids',indices=[100])
-    print '         ++++ Element 200s RIDs ', mat_RID.get_values(location='centroids',indices=[200])
-    print '         ++++ Element 300s RIDs ', mat_RID.get_values(location='centroids',indices=[300])
-    
-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)
+print model_data.basename+' >>>> Creating and assigning surface roughness data for use in evolve'
+
+N = len(domain)
+material_variables = zeros((N, 5), Float)
+
+# Set default material variables n0, d1, n1, d2, n2
+material_variables[:,0] = 0.035 # n0
+material_variables[:,1] = 0.3   # d1  
+material_variables[:,2] = 0.050 # n1  
+material_variables[:,3] = 1.5   # d2  
+material_variables[:,4] = 0.030 # n2  
+
+
+# Assign polygonal data to material_variables from model_data.mat_poly_list
+
+
+points = domain.get_centroid_coordinates(absolute=True)
+
+for k, poly in enumerate(model_data.mat_poly_list):
+    
+    # Get indices of triangles inside polygon k
+    indices = inside_polygon(points, poly, closed=True)
+
+    # Index into material data for polygon k
+    rid = model_data.mat_RID_list[k]
+        
+    # Material data for polygon     
+    n0, d1, n1, d2, n2 = model_data.mat_roughness_data[rid]    
+        
+    for i in indices:
+        material_variables[i, 0] = n0
+        material_variables[i, 1] = d1        
+        material_variables[i, 2] = n1
+        material_variables[i, 3] = d2
+        material_variables[i, 4] = n2                        
+        
+    
+
+print '         ++++ Upgraded material_variables spatially with the material values set in the mif/mid files ' 
+
+    
+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)
 
 print model_data.basename+' >>>> Completed domain construction at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
+
 
 
@@ -316 +329 @@
     print '         ++++ Recomputing friction at t= ', t
         
-    noelementsupdated = 0                                                      # intialise update check
+
     # rebuild the 'friction' values adjusted for depth 
     for i in domain.get_wet_elements():                                        # loop for each wet element in domain
     
-        element_RID_value=mat_RID.get_values(location='centroids',
-                                             indices=[i]) # get current elements centroidal RID value
-                                             
-        # print '         ++++ Element %i centroidal element_RID_value is %s' %( i,element_RID_value[0] )
-        for k in range(len(model_data.mat_roughness_data_list)):               # loop for each matl in tmf list 
-            if int(element_RID_value[0]) == int( model_data.mat_roughness_data_list[k][0] )  :  # got tmf line for that matl
-    #            print '         ++++ found tmf match for RID %i of %s ' %( int(element_RID_value[0]),model_data.mat_roughness_data_list[k][8].rstrip(' ')  )
-                d1 = float(model_data.mat_roughness_data_list[k][4])
-                n1 = float(model_data.mat_roughness_data_list[k][5])
-                d2 = float(model_data.mat_roughness_data_list[k][6])
-                n2 = float(model_data.mat_roughness_data_list[k][7])
-                noelementsupdated = noelementsupdated + 1                       # increment update counter start at zeroth element
-                if   depth.get_values(location='centroids',
-                                      indices=[i]) <= d1 : # recompute friction values from depth for this element get_tuflow_data.py
-
-                    domain.set_quantity('friction',
-                                        numeric=n1,
-                                        location='centroids',
-                                        indices=[i],
-                                        verbose=model_data.anuga_verbose)
-                                        
-                elif depth.get_values(location='centroids',
-                                      indices=[i]) >= d2 :
-                    domain.set_quantity('friction',numeric=n2,location='centroids',
-                                        indices=[i],verbose=model_data.anuga_verbose)
-                else :
-                    n3 = n1+(n2-n1)/(d2-d1)* ( depth.get_values(location='centroids',indices=[i])[0] )
-                    domain.set_quantity('friction',numeric=n3,location='centroids',
-                                        indices=[i],verbose=model_data.anuga_verbose)
-
+        # Get roughness variables 
+        n0 = material_variables[i,0]
+        d1 = material_variables[i,1]
+        n1 = material_variables[i,2]
+        d2 = material_variables[i,3]
+        n2 = material_variables[i,4]                
+                
+
+        # Recompute friction values from depth for this element get_tuflow_data.py
+        d = depth.get_values(location='centroids', indices=[i])[0]
+         
+        if d <= d1: 
+            depth_dependent_friction = n1
+        elif d >= d2:
+            depth_dependent_friction = n2                        
+        else:
+            depth_dependent_friction = n1+(n2-n1)/(d2-d1)*d                         
+
+
+        domain.set_quantity('friction',
+                            numeric=depth_dependent_friction,
+                            location='centroids',
+                            indices=[i],
+                            verbose=model_data.anuga_verbose)                                
+              
               
     if model_data.model_verbose :
-    
         friction = domain.get_quantity('friction')
         
-        print "         ++++  Check  some  RID,depth, friction values returned during this timestep "   
-        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])
-        
-        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])
-        
-        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] )
-                        
- 
-
+        # Print something
+