Changeset 9596
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/anuga/structures/riverwall.py
r9501 r9596 578 578 """ 579 579 580 domain =self.domain580 domain = self.domain 581 581 582 582 # Preliminary definitions 583 isConnected =True584 printInfo =''583 isConnected = True 584 printInfo = '' 585 585 586 586 if(len(self.names)==0): 587 587 if(verbose): 588 printInfo =printInfo+' There are no riverwalls (P'+str(myid)+')\n'588 printInfo = printInfo+' There are no riverwalls (P'+str(myid)+')\n' 589 589 return [printInfo, isConnected] 590 590 591 591 # Shorthand notation 592 rwd =self592 rwd = self 593 593 594 594 for i, name in enumerate(rwd.names): 595 595 # Get indices of edges on this riverwall 596 riverwalledgeInds =rwd.riverwall_edges[(rwd.hydraulic_properties_rowIndex==i).nonzero()[0]]596 riverwalledgeInds = rwd.riverwall_edges[(rwd.hydraulic_properties_rowIndex==i).nonzero()[0]] 597 597 598 598 if(len(riverwalledgeInds)==0): 599 printInfo =printInfo+'Riverwall '+name+' was not found on this mesh (if this is wrong, adjust tol in create_riverwalls)\n'599 printInfo = printInfo+'Riverwall '+name+' was not found on this mesh (if this is wrong, adjust tol in create_riverwalls)\n' 600 600 continue 601 601 # Get their corresponding vertices … … 604 604 # Flag telling us if vertex points are on the boundary of the model 605 605 # Used to help detect disconnected riverwalls (due to round-off) 606 v1_on_boundary =rwd.is_vertex_on_boundary(riverwallV1Inds)607 v2_on_boundary =rwd.is_vertex_on_boundary(riverwallV2Inds)606 v1_on_boundary = rwd.is_vertex_on_boundary(riverwallV1Inds) 607 v2_on_boundary = rwd.is_vertex_on_boundary(riverwallV2Inds) 608 608 609 609 # With discontinuous triangles, we expect edges to occur twice 610 610 # Let's remove duplicates to simplify the analysis 611 repeat =riverwalledgeInds*0612 lre =len(riverwalledgeInds)611 repeat = riverwalledgeInds*0 612 lre = len(riverwalledgeInds) 613 613 # edge coordinates as a complex number, for easy equality checking 614 complex_edge_coordinates =domain.edge_coordinates[riverwalledgeInds,0]+\615 1j*domain.edge_coordinates[riverwalledgeInds,1]614 complex_edge_coordinates = domain.edge_coordinates[riverwalledgeInds,0]+\ 615 1j*domain.edge_coordinates[riverwalledgeInds,1] 616 616 for j in range(lre-1): 617 617 # Ignore if already checked … … 619 619 continue 620 620 # Check for a dupulicate 621 dups =(complex_edge_coordinates[(j+1):lre]==complex_edge_coordinates[j]).nonzero()[0]621 dups = (complex_edge_coordinates[(j+1):lre]==complex_edge_coordinates[j]).nonzero()[0] 622 622 if(len(dups)>0): 623 repeat[dups+j+1] =1623 repeat[dups+j+1] = 1 624 624 625 unique_riverwall_edge_indices =(repeat==0).nonzero()[0]625 unique_riverwall_edge_indices = (repeat==0).nonzero()[0] 626 626 627 627 # Finally, get 'unqiue' edges in the riverwall 628 uEdges =riverwalledgeInds[unique_riverwall_edge_indices]629 uV1 =riverwallV1Inds[unique_riverwall_edge_indices]630 uV2 =riverwallV2Inds[unique_riverwall_edge_indices]631 uV1_boundary =v1_on_boundary[unique_riverwall_edge_indices]632 uV2_boundary =v2_on_boundary[unique_riverwall_edge_indices]628 uEdges = riverwalledgeInds[unique_riverwall_edge_indices] 629 uV1 = riverwallV1Inds[unique_riverwall_edge_indices] 630 uV2 = riverwallV2Inds[unique_riverwall_edge_indices] 631 uV1_boundary = v1_on_boundary[unique_riverwall_edge_indices] 632 uV2_boundary = v2_on_boundary[unique_riverwall_edge_indices] 633 633 634 634 # Next, count how many times each vertex value occurs. 635 635 # For a 'connected' riverwall, we only want 2 edges where a vertex occurs only once, 636 636 # unless the vertex is on the boundary of the domain 637 lure =len(uEdges)638 complex_v1_coordinates =domain.vertex_coordinates[uV1,0]+\639 1j*domain.vertex_coordinates[uV1,1]640 complex_v2_coordinates =domain.vertex_coordinates[uV2,0]+\641 1j*domain.vertex_coordinates[uV2,1]642 v1Counter =uEdges*0643 v2Counter =uEdges*0637 lure = len(uEdges) 638 complex_v1_coordinates = domain.vertex_coordinates[uV1,0]+\ 639 1j*domain.vertex_coordinates[uV1,1] 640 complex_v2_coordinates = domain.vertex_coordinates[uV2,0]+\ 641 1j*domain.vertex_coordinates[uV2,1] 642 v1Counter = uEdges*0 643 v2Counter = uEdges*0 644 644 for j in range(lure): 645 v1Counter[j] =(complex_v1_coordinates==complex_v1_coordinates[j]).sum()+\646 (complex_v2_coordinates==complex_v1_coordinates[j]).sum()647 v2Counter[j] =(complex_v1_coordinates==complex_v2_coordinates[j]).sum()+\648 (complex_v2_coordinates==complex_v2_coordinates[j]).sum()649 650 num_disconnected_edges =((v1Counter==1)*(1-uV1_boundary)).sum()+\651 ((v2Counter==1)*(1-uV2_boundary)).sum()645 v1Counter[j] = (complex_v1_coordinates==complex_v1_coordinates[j]).sum()+\ 646 (complex_v2_coordinates==complex_v1_coordinates[j]).sum() 647 v2Counter[j] = (complex_v1_coordinates==complex_v2_coordinates[j]).sum()+\ 648 (complex_v2_coordinates==complex_v2_coordinates[j]).sum() 649 650 num_disconnected_edges = ((v1Counter==1)*(1-uV1_boundary)).sum()+\ 651 ((v2Counter==1)*(1-uV2_boundary)).sum() 652 652 653 653 if(verbose): 654 printInfo =printInfo+ ' On riverwall '+ str(name) +' there are '+ str(num_disconnected_edges)+\654 printInfo = printInfo+ ' On riverwall '+ str(name) +' there are '+ str(num_disconnected_edges)+\ 655 655 ' endpoints inside the domain [ignoring points on the boundary polygon] (P'+str(myid)+')\n' 656 656 657 if(num_disconnected_edges <=2):657 if(num_disconnected_edges <= 2): 658 658 if(verbose): 659 659 pass 660 660 #printInfo=printInfo+ " This is consistent with a continuous wall \n" 661 661 else: 662 isConnected =False663 printInfo =printInfo+' Riverwall ' + name +' appears to be discontinuous. (P'+str(myid)+')\n'+\662 isConnected = False 663 printInfo = printInfo + ' Riverwall ' + name +' appears to be discontinuous. (P'+str(myid)+')\n'+\ 664 664 ' This suggests there is a gap in the wall, which should not occur\n' 665 665 … … 683 683 return 684 684 685 if(myid ==0):685 if(myid == 0): 686 686 # Make output directory 687 687 try: … … 691 691 # Make output files with empty contents 692 692 for i, riverWallFile in enumerate(self.names): 693 newFile=open(output_dir +'/'+os.path.splitext(os.path.basename(riverWallFile))[0]+'.txt','w')693 newFile=open(output_dir + '/' + os.path.splitext(os.path.basename(riverWallFile))[0] + '.txt','w') 694 694 # Write hydraulic variable information 695 695 hydraulicVars=self.hydraulic_properties[i,:] 696 696 newFile.write('## Hydraulic Variable values below ## \n') 697 newFile.write(str(self.hydraulic_variable_names) +'\n')698 newFile.write(str(hydraulicVars) +'\n')697 newFile.write(str(self.hydraulic_variable_names) + '\n') 698 newFile.write(str(hydraulicVars) + '\n') 699 699 newFile.write('\n') 700 700 newFile.write('## xyElevation at edges below. Order may be erratic for parallel runs ##\n') … … 716 716 for i in range(numprocs): 717 717 # Write 1 processor at a time 718 if(myid ==i):718 if(myid == i): 719 719 for j, riverWallname in enumerate(self.names): 720 720 # Get xyz data for riverwall j 721 riverWallInds =(self.hydraulic_properties_rowIndex==j).nonzero()[0].tolist()722 riverWallDomainInds =self.riverwall_edges[riverWallInds].tolist()723 myXCoords =domain.edge_coordinates[riverWallDomainInds,0]+domain.geo_reference.xllcorner724 myYCoords =domain.edge_coordinates[riverWallDomainInds,1]+domain.geo_reference.yllcorner725 myElev =self.riverwall_elevation[riverWallInds]721 riverWallInds = (self.hydraulic_properties_rowIndex==j).nonzero()[0].tolist() 722 riverWallDomainInds = self.riverwall_edges[riverWallInds].tolist() 723 myXCoords = domain.edge_coordinates[riverWallDomainInds,0] + domain.geo_reference.xllcorner 724 myYCoords = domain.edge_coordinates[riverWallDomainInds,1] + domain.geo_reference.yllcorner 725 myElev = self.riverwall_elevation[riverWallInds] 726 726 727 727 # Open file for appending data 728 theFile =open(output_dir+'/'+os.path.splitext(os.path.basename(riverWallname))[0]+'.txt','a')728 theFile = open(output_dir + '/' + os.path.splitext(os.path.basename(riverWallname))[0] + '.txt','a') 729 729 for k in range(len(myElev)): 730 theFile.write(str(myXCoords[k]) +','+str(myYCoords[k])+','+str(myElev[k])+'\n')730 theFile.write(str(myXCoords[k]) + ',' + str(myYCoords[k]) + ',' + str(myElev[k]) + '\n') 731 731 theFile.close() 732 732 -
trunk/anuga_core/anuga/structures/structure_operator.py
r9130 r9596 36 36 manning=None, 37 37 enquiry_gap=None, 38 use_momentum_jet=False, 39 zero_outflow_momentum=True, 38 40 description=None, 39 41 label=None, … … 84 86 self.manning = manning 85 87 self.enquiry_gap = enquiry_gap 88 self.use_momentum_jet = use_momentum_jet 89 self.zero_outflow_momentum = zero_outflow_momentum 90 91 if use_momentum_jet and zero_outflow_momentum: 92 msg = "Can't have use_momentum_jet and zero_outflow_momentum both True" 93 raise Exception(msg) 86 94 87 95 if description == None: … … 222 230 223 231 loss = (old_inflow_depth - new_inflow_depth)*self.inflow.get_area() 232 xmom_loss = (old_inflow_xmom - new_inflow_xmom)*self.inflow.get_area() 233 ymom_loss = (old_inflow_ymom - new_inflow_ymom)*self.inflow.get_area() 224 234 225 235 # set outflow … … 229 239 timestep_star = 0.0 230 240 231 232 233 234 241 outflow_extra_depth = Q*timestep_star/self.outflow.get_area() 235 242 outflow_direction = - self.outflow.outward_culvert_vector 236 outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction243 #outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction 237 244 238 245 gain = outflow_extra_depth*self.outflow.get_area() … … 262 269 outflow_extra_depth = Q*timestep/self.outflow.get_area() 263 270 outflow_direction = - self.outflow.outward_culvert_vector 264 outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction271 #outflow_extra_momentum = outflow_extra_depth*barrel_speed*outflow_direction 265 272 266 273 gain = outflow_extra_depth*self.outflow.get_area() … … 275 282 276 283 self.accumulated_flow += gain 277 self.discharge = outflow_extra_depth*self.outflow.get_area()/timestep #Q284 self.discharge = Q*timestep_star/timestep # outflow_extra_depth*self.outflow.get_area()/timestep #Q 278 285 self.velocity = barrel_speed # self.discharge/outlet_depth/self.width 279 286 self.outlet_depth = outlet_depth … … 281 288 new_outflow_depth = self.outflow.get_average_depth() + outflow_extra_depth 282 289 283 if self.use_momentum_jet : 290 self.outflow.set_depths(new_outflow_depth) 291 292 if self.use_momentum_jet: 284 293 # FIXME (SR) Review momentum to account for possible hydraulic jumps at outlet 294 # FIXME (GD) Depending on barrel speed I think this will be either 295 # a source or sink of momentum (considering the momentum losses 296 # above). Might not always be reasonable. 285 297 #new_outflow_xmom = self.outflow.get_average_xmom() + outflow_extra_momentum[0] 286 298 #new_outflow_ymom = self.outflow.get_average_ymom() + outflow_extra_momentum[1] 287 299 new_outflow_xmom = barrel_speed*new_outflow_depth*outflow_direction[0] 288 300 new_outflow_ymom = barrel_speed*new_outflow_depth*outflow_direction[1] 289 290 else: 301 302 elif self.zero_outflow_momentum: 303 # Default case (FIXME: Should this be default? We lose all the 304 # momentum, often not physically reasonable) 305 new_outflow_xmom = 0.0 306 new_outflow_ymom = 0.0 291 307 #new_outflow_xmom = outflow.get_average_xmom() 292 308 #new_outflow_ymom = outflow.get_average_ymom() 293 309 294 new_outflow_xmom = 0.0 295 new_outflow_ymom = 0.0 296 297 self.outflow.set_depths(new_outflow_depth) 310 else: 311 # Add the momentum lost from the inflow to the outflow. For 312 # structures where barrel_speed is unknown + direction doesn't 313 # change from inflow to outflow 314 new_outflow_xmom = self.outflow.get_average_xmom() + xmom_loss/self.outflow.get_area() 315 new_outflow_ymom = self.outflow.get_average_ymom() + ymom_loss/self.outflow.get_area() 316 298 317 self.outflow.set_xmoms(new_outflow_xmom) 299 318 self.outflow.set_ymoms(new_outflow_ymom) … … 514 533 self.log_filename = self.label + '.log' 515 534 log_to_file(self.log_filename, self.statistics(), mode='w') 516 log_to_file(self.log_filename, 'time, discharge,velocity,driving_energy,delta_total_energy')535 log_to_file(self.log_filename, 'time, discharge, velocity, accumulated_flow, driving_energy, delta_total_energy') 517 536 518 537 #log_to_file(self.log_filename, self.culvert_type)
Note: See TracChangeset
for help on using the changeset viewer.
