Changeset 7971
 Timestamp:
 Aug 24, 2010, 11:06:01 PM (12 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7968 r7971 6 6 from anuga.geometry.polygon import plot_polygons, polygon_area 7 7 8 9 8 from anuga.utilities.numerical_tools import mean 10 9 from anuga.utilities.numerical_tools import ensure_numeric, sign … … 14 13 import anuga.utilities.log as log 15 14 15 import Inlet 16 16 17 import numpy as num 17 from math import sqrt 18 import math 18 19 19 20 class Below_interval(Exception): pass 20 21 class Above_interval(Exception): pass 21 22 23 24 class Inlet:25 """Contains information associated with each inlet26 """27 28 def __init__(self,domain,polygon,enquiry_point,inlet_vector):29 30 self.domain = domain31 self.domain_bounding_polygon = self.domain.get_boundary_polygon()32 self.polygon = polygon33 self.enquiry_point = enquiry_point34 self.inlet_vector = inlet_vector35 36 # FIXME (SR) Using get_triangle_containing_point which needs to be sped up37 self.enquiry_index = self.domain.get_triangle_containing_point(self.enquiry_point)38 39 self.compute_inlet_triangle_indices()40 41 42 def compute_inlet_averages(self):43 44 45 46 self.cell_indices = self.exchange_triangle_indices[0]47 self.areas = areas[cell_indices]48 49 50 # Inlet Averages51 self.heights = stage[cell_indices]elevation[cell_indices]52 self.total_water_volume = num.sum(self.heights*self.areas)53 self.average_height = self.total_water_volume/self.total_inlet_area54 55 56 57 58 def compute_inlet_triangle_indices(self):59 60 # Get boundary (in absolute coordinates)61 bounding_polygon = self.domain_bounding_polygon62 centroids = self.domain.get_centroid_coordinates(absolute=True)63 64 inlet_polygon = self.polygon65 66 # Check that polygon lies within the mesh.67 for point in inlet_polygon:68 msg = 'Point %s in polygon for forcing term' % str(point)69 msg += ' did not fall within the domain boundary.'70 assert is_inside_polygon(point, bounding_polygon), msg71 72 73 self.triangle_indices = inside_polygon(centroids, inlet_polygon)74 75 if len(self.triangle_indices) == 0:76 region = 'Inlet polygon=%s' % (inlet_polygon)77 msg = 'No triangles have been identified in '78 msg += 'specified region: %s' % region79 raise Exception, msg80 81 # Compute exchange area as the sum of areas of triangles identified82 # by polygon83 self.area = 0.084 for j in self.triangle_indices:85 self.area += self.domain.areas[j]86 87 88 msg = 'Inlet exchange area has area = %f' % self.area89 assert self.area > 0.090 91 22 92 23 … … 117 48 self.domain.set_fractional_step_operator(self) 118 49 119 self.end_points = [end_point0, end_point1]50 self.end_points = [end_point0, end_point1] 120 51 self.enquiry_gap_factor = enquiry_gap_factor 121 52 … … 134 65 #FIXME (SR) Put this into a foe loop to deal with more inlets 135 66 self.inlets = [] 136 polygon0 = self. exchange_polygons[0]67 polygon0 = self.inlet_polygons[0] 137 68 enquiry_pt0 = self.enquiry_points[0] 138 69 inlet0_vector = self.culvert_vector 139 140 self.inlets.append(Inlet(self.domain,polygon0,enquiry_pt0,inlet0_vector)) 141 142 polygon1 = self.exchange_polygons[1] 70 self.inlets.append(Inlet.Inlet(self.domain, polygon0, enquiry_pt0, inlet0_vector)) 71 72 polygon1 = self.inlet_polygons[1] 143 73 enquiry_pt1 = self.enquiry_points[1] 144 74 inlet1_vector =  self.culvert_vector 145 146 self.inlets.append(Inlet(self.domain,polygon1,enquiry_pt1, inlet1_vector)) 147 148 149 # aliases to quantity centroid values and cell areas 150 self.areas = self.domain.areas 151 self.stage = self.domain.quantities['stage'].centroid_values 152 self.elevation = self.domain.quantities['elevation'].centroid_values 153 self.xmom = self.domain.quantities['xmomentum'].centroid_values 154 self.ymom = self.domain.quantities['ymomentum'].centroid_values 155 75 self.inlets.append(Inlet.Inlet(self.domain, polygon1, enquiry_pt1, inlet1_vector)) 156 76 77 self.areas = self.domain.areas 78 self.stages = self.domain.quantities['stage'].centroid_values 79 self.elevations = self.domain.quantities['elevation'].centroid_values 80 self.xmoms = self.domain.quantities['xmomentum'].centroid_values 81 self.ymoms = self.domain.quantities['ymomentum'].centroid_values 82 157 83 self.print_stats() 158 84 159 85 160 161 86 def __call__(self): 162 163 87 164 88 # Time stuff … … 166 90 timestep = self.domain.get_timestep() 167 91 168 169 92 inlet0 = self.inlets[0] 170 93 inlet1 = self.inlets[1] 171 172 94 173 95 # Aliases to cell indices … … 175 97 inlet1_indices = inlet1.triangle_indices 176 98 177 178 99 # Inlet0 averages 179 inlet0_heights = self.stage[inlet0_indices]self.elevation[inlet0_indices] 180 inlet0_areas = self.areas[inlet0_indices] 181 100 inlet0_heights = inlet0.get_heights() 101 inlet0_areas = inlet0.get_areas() 182 102 inlet0_water_volume = num.sum(inlet0_heights*inlet0_areas) 183 184 103 average_inlet0_height = inlet0_water_volume/inlet0.area 185 104 186 105 # Inlet1 averages 187 inlet1_heights = self.stage[inlet1_indices]self.elevation[inlet1_indices] 188 inlet1_areas = self.areas[inlet1_indices] 189 106 inlet1_heights = inlet1.get_heights() 107 inlet1_areas = inlet1.get_areas() 190 108 inlet1_water_volume = num.sum(inlet1_heights*inlet1_areas) 191 192 109 average_inlet1_height = inlet1_water_volume/inlet1.area 193 194 110 195 111 # Transfer 196 112 transfer_water = timestep*inlet0_water_volume 197 198 self.stage[inlet0_indices] = self.elevation[inlet0_indices] + average_inlet0_height  transfer_water 199 self.xmom[inlet0_indices] = 0.0 200 self.ymom[inlet0_indices] = 0.0 201 202 203 self.stage[inlet1_indices] = self.elevation[inlet1_indices] + average_inlet1_height + transfer_water 204 self.xmom[inlet1_indices] = 0.0 205 self.ymom[inlet1_indices] = 0.0 113 114 self.stages.put(inlet0_indices, inlet0.get_elevations() + average_inlet0_height  transfer_water) 115 self.xmoms.put(inlet0_indices, 0.0) 116 self.ymoms.put(inlet0_indices, 0.0) 117 118 self.stages.put(inlet1_indices, inlet1.get_elevations() + average_inlet1_height + transfer_water) 119 self.xmoms.put(inlet1_indices, 0.0) 120 self.ymoms.put(inlet1_indices, 0.0) 206 121 207 122 … … 232 147 233 148 234 235 236 237 238 239 149 def set_store_hydrograph_discharge(self, filename=None): 240 150 … … 251 161 fid.close() 252 162 163 253 164 def create_culvert_polygons(self): 165 254 166 """Create polygons at the end of a culvert inlet and outlet. 255 167 At either end two polygons will be created; one for the actual flow to pass through and one a little further away … … 265 177 266 178 self.culvert_vector = num.array([dx, dy]) 267 self.culvert_length = sqrt(num.sum(self.culvert_vector**2))179 self.culvert_length = math.sqrt(num.sum(self.culvert_vector**2)) 268 180 assert self.culvert_length > 0.0, 'The length of culvert is less than 0' 269 181 … … 278 190 gap = (1 + self.enquiry_gap_factor)*h 279 191 280 self. exchange_polygons = []192 self.inlet_polygons = [] 281 193 self.enquiry_points = [] 282 194 … … 288 200 p2 = p1 + i0*h 289 201 p3 = p0 + i0*h 290 self. exchange_polygons.append(num.array([p0, p1, p2, p3]))202 self.inlet_polygons.append(num.array([p0, p1, p2, p3])) 291 203 self.enquiry_points.append(self.end_points[i] + i0*gap) 292 204 293 294 295 296 # Check that enquiry points are outside exchange polygons 205 # Check that enquiry points are outside inlet polygons 297 206 for i in [0,1]: 298 polygon = self. exchange_polygons[i]207 polygon = self.inlet_polygons[i] 299 208 # FIXME (SR) Probably should calculate the area of all the triangles 300 209 # associated with this polygon, as there is likely to be some … … 313 222 assert not inside_polygon(point, polygon), msg 314 223 315 316 317 318 319 320 321 322 224 323 225 def adjust_flow_for_available_water_at_inlet(self, Q, delta_t): 324 226 """Adjust Q downwards depending on available water at inlet … … 463 365 """Compute new rates for inlet and outlet 464 366 """ 465 466 # Short hands467 domain = self.domain468 dq = domain.quantities469 470 367 # Time stuff 471 time = domain.get_time()368 time = self.domain.get_time() 472 369 self.last_update = time 473 370 474 475 371 if hasattr(self, 'log_filename'): 476 372 log_filename = self.log_filename … … 478 374 # Compute stage, energy and velocity at the 479 375 # enquiry points at each end of the culvert 480 openings = self.openings 481 for i, opening in enumerate(openings): 482 idx = self.enquiry_indices[i] 483 484 stage = dq['stage'].get_values(location='centroids', 485 indices=[idx])[0] 486 depth = h = stageopening.elevation 487 488 489 # Get velocity 490 xmomentum = dq['xmomentum'].get_values(location='centroids', 491 indices=[idx])[0] 492 ymomentum = dq['xmomentum'].get_values(location='centroids', 493 indices=[idx])[0] 494 495 if h > minimum_allowed_height: 496 u = xmomentum/(h + velocity_protection/h) 497 v = ymomentum/(h + velocity_protection/h) 376 total_energies = [] 377 specific_energies = [] 378 velocities = [] 379 380 for i, inlet in enumerate(self.inlets): 381 382 stage = inlet.get_average_stage() 383 depth = stage  inlet.get_average_elevation() 384 385 if depth > minimum_allowed_height: 386 u, v = inlet.get_average_velocities 498 387 else: 499 u = v = 0.0 388 u = 0.0 389 v = 0.0 500 390 501 v_squared = u*u + v*v502 391 v_squared = u**2 + v**2 392 503 393 if self.use_velocity_head is True: 504 394 velocity_head = 0.5*v_squared/g … … 506 396 velocity_head = 0.0 507 397 508 opening.total_energy = velocity_head + stage 509 opening.specific_energy = velocity_head + depth 510 opening.stage = stage 511 opening.depth = depth 512 opening.velocity = sqrt(v_squared) 513 398 total_energies[i] = velocity_head + stage 399 specific_energies[i] = velocity_head + depth 400 velocities[i] = math.sqrt(v_squared) 514 401 515 402 # We now need to deal with each opening individually 516 403 # Determine flow direction based on total energy difference 517 delta_total_energy = openings[0].total_energy  openings[1].total_energy 518 if delta_total_energy > 0: 519 inlet = openings[0] 520 outlet = openings[1] 521 522 # FIXME: I think this whole momentum jet thing could be a bit more elegant 523 inlet.momentum = self.opening_momentum[0] 524 outlet.momentum = self.opening_momentum[1] 404 delta_total_energy = total_energies[0]  total_energies[1] 405 if delta_total_energy >= 0: 406 inlet = inlets[0] 407 outlet = inlets[1] 525 408 else: 526 inlet = openings[1] 527 outlet = openings[0] 528 529 inlet.momentum = self.opening_momentum[1] 530 outlet.momentum = self.opening_momentum[0] 531 532 delta_total_energy = delta_total_energy 533 534 self.inlet = inlet 535 self.outlet = outlet 536 537 msg = 'Total energy difference is negative' 538 assert delta_total_energy >= 0.0, msg 409 msg = 'Total energy difference is negative' 410 assert delta_total_energy >= 0.0, msg 539 411 540 412 # Recompute slope and issue warning if flow is uphill 541 413 # These values do not enter the computation 542 delta_z = inlet. elevation  outlet.elevation543 culvert_slope = (delta_z/self. length)414 delta_z = inlet.get_average_elevation()  outlet.get_average_elevation() 415 culvert_slope = (delta_z/self.culvert_length) 544 416 if culvert_slope < 0.0: 545 417 # Adverse gradient  flow is running uphill … … 554 426 s = 'Delta total energy = %.3f' %(delta_total_energy) 555 427 log_to_file(log_filename, s) 556 557 428 558 429 # Determine controlling energy (driving head) for culvert … … 564 435 driving_head = inlet.specific_energy 565 436 566 567 568 437 if self.inlet.depth <= self.trigger_depth: 569 438 Q = 0.0 … … 617 486 log_filename=self.log_filename) 618 487 619 620 621 488 # Adjust discharge for multiple barrels 622 489 Q *= self.number_of_barrels … … 628 495 self.outlet.rate = Q 629 496 630 631 497 # Momentum jet stuff 632 498 if self.use_momentum_jet is True: 633 634 499 635 500 # Compute barrel momentum … … 646 511 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y) 647 512 log_to_file(self.log_filename, s) 648 649 513 650 514 # Update momentum
Note: See TracChangeset
for help on using the changeset viewer.