Changeset 5975


Ignore:
Timestamp:
Nov 19, 2008, 1:28:49 PM (15 years ago)
Author:
ole
Message:

Restricted main loop to wet elements

File:
1 edited

Legend:

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

    r5974 r5975  
    7676
    7777
    78 
     78domain
     79    domain.forcing_terms.append(hydrograph)
     80
     81print model_data.basename+' >>>> Completed assignment of %i inflow  hydrographs at t= %.2f hours ' %(len(model_data.infl
    7982print model_data.basename+' >>>> Completed mesh creation at t= %.2f hours ' %(float(time.time()-t0)/3600.0)
    8083
     
    306309for t in domain.evolve(yieldstep = steptime, finaltime = endtime):
    307310   
    308     domain.write_time()                                                            # confirm t each steptime
    309     domain.write_boundary_statistics(tags='Lake', quantities='stage')              # indicate bdry stats each steptime
    310     depth = domain.create_quantity_from_expression('stage - elevation')            # create depth instance for this timestep
    311 
    312     if int(float((t+0.00001)/steptime)) > lastprintstep  :                                   # Recompute friction at each print timestep
    313         lastprintstep = int(float((t+0.00001)/steptime))
    314         print '         ++++ Recomputing friction at t= ', t
    315        
    316         noelementsupdated = 0                                                      # intialise update check
    317         # rebuild the 'friction' values adjusted for depth
    318         for i in range(len(domain)) :                                              # loop for each element in domain
    319        
    320             element_RID_value=mat_RID.get_values(location='centroids',indices=[i]) # get current elements centroidal RID value
    321         #    print '         ++++ Element %i centroidal element_RID_value is %s' %( i,element_RID_value[0] )
    322             for k in range(len(model_data.mat_roughness_data_list)):               # loop for each matl in tmf list
    323                 if int(element_RID_value[0]) == int( model_data.mat_roughness_data_list[k][0] )  :  # got tmf line for that matl
    324         #            print '         ++++ found tmf match for RID %i of %s ' %( int(element_RID_value[0]),model_data.mat_roughness_data_list[k][8].rstrip(' ')  )
    325                     d1 = float(model_data.mat_roughness_data_list[k][4])
    326                     n1 = float(model_data.mat_roughness_data_list[k][5])
    327                     d2 = float(model_data.mat_roughness_data_list[k][6])
    328                     n2 = float(model_data.mat_roughness_data_list[k][7])
    329                     noelementsupdated = noelementsupdated +1                        # increment update counter start at zeroth element
    330                     if   depth.get_values(location='centroids',indices=[i]) <= d1 : # recompute friction values from depth for this element get_tuflow_data.py
    331 
    332                         domain.set_quantity('friction',numeric=n1,location='centroids',indices=[i],verbose=model_data.anuga_verbose)
    333                     elif depth.get_values(location='centroids',indices=[i]) >= d2 :
    334                         domain.set_quantity('friction',numeric=n2,location='centroids',indices=[i],verbose=model_data.anuga_verbose)
    335                     else :
    336                         n3 = n1+(n2-n1)/(d2-d1)* ( depth.get_values(location='centroids',indices=[i])[0] )
    337                         domain.set_quantity('friction',numeric=n3,location='centroids',indices=[i],verbose=model_data.anuga_verbose)
    338 
    339         if i <> (noelementsupdated-1) :                                              # Check all elements actually updated
    340             # EHR - will need to modify if only wet elements updated
    341             print '         ++++ WARNING - %i elements in domain but only %i updated ' %(i+1,noelementsupdated)
     311    domain.write_time()                                                  # confirm t each steptime
     312    domain.write_boundary_statistics(tags='Lake', quantities='stage')    # indicate bdry stats each steptime
     313    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
     314
     315
     316    print '         ++++ Recomputing friction at t= ', t
     317       
     318    noelementsupdated = 0                                                      # intialise update check
     319    # rebuild the 'friction' values adjusted for depth
     320    for i in domain.get_wet_elements():                                        # loop for each wet element in domain
     321   
     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
    342352             
    343353    if model_data.model_verbose :
     354   
     355        friction = domain.get_quantity('friction')
     356       
    344357        print "         ++++  Check  some  RID,depth, friction values returned during this timestep "   
    345         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] )
    346         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] )
    347         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] )
     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] )
    348363                       
    349364 
Note: See TracChangeset for help on using the changeset viewer.