Changeset 5975

Timestamp:

Nov 19, 2008, 1:28:49 PM (16 years ago)

Author:

ole

Message:

Restricted main loop to wet elements

File:

• anuga_work/development/ted_rigby/model_run.py (modified) (2 diffs)

anuga_work/development/ted_rigby/model_run.py

@@ -76 +76 @@
 
 
-
+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)
 
@@ -306 +309 @@
 for t in domain.evolve(yieldstep = steptime, finaltime = endtime):
     
-    domain.write_time()                                                            # confirm t each steptime
-    domain.write_boundary_statistics(tags='Lake', quantities='stage')              # indicate bdry stats each steptime
-    depth = domain.create_quantity_from_expression('stage - elevation')            # create depth instance for this timestep
-
-    if int(float((t+0.00001)/steptime)) > lastprintstep  :                                   # Recompute friction at each print timestep
-        lastprintstep = int(float((t+0.00001)/steptime))
-        print '         ++++ Recomputing friction at t= ', t
-        
-        noelementsupdated = 0                                                      # intialise update check
-        # rebuild the 'friction' values adjusted for depth 
-        for i in range(len(domain)) :                                              # loop for each 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)
-
-        if i <> (noelementsupdated-1) :                                              # Check all elements actually updated
-            # EHR - will need to modify if only wet elements updated
-            print '         ++++ WARNING - %i elements in domain but only %i updated ' %(i+1,noelementsupdated)
+    domain.write_time()                                                  # confirm t each steptime
+    domain.write_boundary_statistics(tags='Lake', quantities='stage')    # indicate bdry stats each steptime
+    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
+
+
+    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)
+
               
     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]),domain.get_quantity('friction',location='centroids',indices=[100] )
-        print "               element 200 ",mat_RID.get_values(location='centroids',indices=[200]),depth.get_values(location='centroids',indices=[200]),domain.get_quantity('friction',location='centroids',indices=[200] )
-        print "               element 300 ",mat_RID.get_values(location='centroids',indices=[300]),depth.get_values(location='centroids',indices=[300]),domain.get_quantity('friction',location='centroids',indices=[300] )
+        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] )