Changeset 7975
- Timestamp:
- Aug 25, 2010, 4:35:40 PM (14 years ago)
- Location:
- trunk/anuga_core/source/anuga/structures
- Files:
-
- 3 edited
- 1 moved
-
boyd_box_routine.py (modified) (1 diff)
-
culvert_operator.py (modified) (5 diffs)
-
inlet.py (moved) (moved from trunk/anuga_core/source/anuga/structures/Inlet.py) (4 diffs)
-
testing_culvert.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/anuga_core/source/anuga/structures/boyd_box_routine.py
r7969 r7975 9 9 10 10 11 def boyd_box( height, width, flow_width, inflow_specific_energy):11 def boyd_box(in): 12 12 """Boyd's generalisation of the US department of transportation culvert methods 13 13 -
trunk/anuga_core/source/anuga/structures/culvert_operator.py
r7974 r7975 13 13 import anuga.utilities.log as log 14 14 15 import Inlet15 import inlet 16 16 17 17 import numpy as num … … 68 68 enquiry_pt0 = self.enquiry_points[0] 69 69 inlet0_vector = self.culvert_vector 70 self.inlets.append( Inlet.Inlet(self.domain, polygon0, enquiry_pt0, inlet0_vector))70 self.inlets.append(inlet.Inlet(self.domain, polygon0, enquiry_pt0, inlet0_vector)) 71 71 72 72 polygon1 = self.inlet_polygons[1] 73 73 enquiry_pt1 = self.enquiry_points[1] 74 74 inlet1_vector = - self.culvert_vector 75 self.inlets.append( Inlet.Inlet(self.domain, polygon1, enquiry_pt1, inlet1_vector))75 self.inlets.append(inlet.Inlet(self.domain, polygon1, enquiry_pt1, inlet1_vector)) 76 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 77 82 78 83 79 self.print_stats() … … 90 86 timestep = self.domain.get_timestep() 91 87 92 inlet0 = self.inlets[0] 93 inlet1 = self.inlets[1] 94 95 # Aliases to cell indices 96 inlet0_indices = inlet0.triangle_indices 97 inlet1_indices = inlet1.triangle_indices 98 99 # Inlet0 averages 100 inlet0_heights = inlet0.get_heights() 101 inlet0_areas = inlet0.get_areas() 102 inlet0_water_volume = num.sum(inlet0_heights*inlet0_areas) 103 average_inlet0_height = inlet0_water_volume/inlet0.area 104 105 # Inlet1 averages 106 inlet1_heights = inlet1.get_heights() 107 inlet1_areas = inlet1.get_areas() 108 inlet1_water_volume = num.sum(inlet1_heights*inlet1_areas) 109 average_inlet1_height = inlet1_water_volume/inlet1.area 88 89 90 inflow = self.inlets[0] 91 outflow = self.inlets[1] 92 93 # Determine flow direction based on total energy difference 94 delta_total_energy = inflow.get_average_total_energy() - outflow.get_average_total_energy() 95 96 if delta_total_energy < 0: 97 inflow = self.inlets[1] 98 outflow = self.inlets[0] 99 delta_total_energy = -delta_total_energy 100 101 delta_z = inflow.get_average_elevation() - outflow.get_average_elevation() 102 culvert_slope = delta_z/self.culvert_length 103 104 # Determine controlling energy (driving head) for culvert 105 if inflow.get_average_specific_energy() > delta_total_energy: 106 # Outlet control 107 driving_head = delta_total_energy 108 else: 109 # Inlet control 110 driving_head = inflow.get_average_specific_energy() 111 110 112 111 113 # Transfer 112 transfer_water = timestep*inlet0_water_volume 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) 114 from culvert_routines import boyd_generalised_culvert_model 115 Q, barrel_velocity, culvert_outlet_depth =\ 116 boyd_generalised_culvert_model(inflow.get_average_height(), 117 outflow.get_average_height(), 118 inflow.get_average_speed(), 119 outflow.get_average_speed(), 120 inflow.get_average_specific_energy(), 121 delta_total_energy, 122 g, 123 culvert_length=self.culvert_length, 124 culvert_width=self.width, 125 culvert_height=self.height, 126 culvert_type='box', 127 manning=0.01) 128 129 transfer_water = Q*timestep 130 131 132 inflow.set_heights(inflow.get_average_height() - transfer_water) 133 inflow.set_xmoms(0.0) 134 inflow.set_ymoms(0.0) 135 136 137 outflow.set_heights(outflow.get_average_height() + transfer_water) 138 outflow.set_xmoms(0.0) 139 outflow.set_ymoms(0.0) 140 121 141 122 142 … … 147 167 148 168 149 def set_store_hydrograph_discharge(self, filename=None): 150 151 if filename is None: 152 self.filename = 'culvert_discharge_hydrograph' 153 else: 154 self.filename = filename 155 156 self.discharge_hydrograph = True 157 158 self.timeseries_filename = self.filename + '_timeseries.csv' 159 fid = open(self.timeseries_filename, 'w') 160 fid.write('time, discharge\n') 161 fid.close() 169 162 170 163 171 … … 223 231 224 232 225 def adjust_flow_for_available_water_at_inlet(self, Q, delta_t): 226 """Adjust Q downwards depending on available water at inlet 227 228 This is a critical step in modelling bridges and Culverts 229 the predicted flow through a structure based on an abstract 230 algorithm can at times request for water that is simply not 231 available due to any number of constrictions that limit the 232 flow approaching the structure In order to ensure that 233 there is adequate flow available certain checks are 234 required There needs to be a check using the Static Water 235 Volume sitting infront of the structure, In addition if the 236 water is moving the available water will be larger than the 237 static volume 238 239 NOTE To temporarily switch this off for Debugging purposes 240 rem out line in function def compute_rates below 241 """ 242 243 if delta_t < epsilon: 244 # No need to adjust if time step is very small or zero 245 # In this case the possible flow will be very large 246 # anyway. 247 return Q 248 249 # Short hands 250 domain = self.domain 251 dq = domain.quantities 252 time = domain.get_time() 253 I = self.inlet 254 idx = I.exchange_indices 255 256 # Find triangle with the smallest depth 257 stage = dq['stage'].get_values(location='centroids', 258 indices=[idx]) 259 elevation = dq['elevation'].get_values(location='centroids', 260 indices=[idx]) 261 depth = stage-elevation 262 min_depth = min(depth.flat) # This may lead to errors if edge of area is at a higher level !!!! 263 avg_depth = mean(depth.flat) # Yes, but this one violates the conservation unit tests 264 265 266 267 # FIXME (Ole): If you want these, use log.critical() and 268 # make the statements depend on verbose 269 #print I.depth 270 #print I.velocity 271 #print self.width 272 273 # max_Q Based on Volume Calcs 274 275 276 depth_term = min_depth*I.exchange_area/delta_t 277 if min_depth < 0.2: 278 # Only add velocity term in shallow waters (< 20 cm) 279 # This is a little ad hoc, but maybe it is reasonable 280 velocity_term = self.width*min_depth*I.velocity 281 else: 282 velocity_term = 0.0 283 284 # This one takes approaching water into account 285 max_Q = max(velocity_term, depth_term) 286 287 # This one preserves Volume 288 #max_Q = depth_term 289 290 291 if self.verbose is True: 292 log.critical('Max_Q = %f' % max_Q) 293 msg = 'Width = %.2fm, Depth at inlet = %.2f m, Velocity = %.2f m/s. ' % (self.width, I.depth, I.velocity) 294 msg += 'Max Q = %.2f m^3/s' %(max_Q) 295 log.critical(msg) 296 297 if self.log_filename is not None: 298 log_to_file(self.log_filename, msg) 299 # New Procedure for assessing the flow available to the Culvert 300 # This routine uses the GET FLOW THROUGH CROSS SECTION 301 # Need to check Several Polyline however 302 # Firstly 3 sides of the exchange Poly 303 # then only the Line Directly infront of the Polygon 304 # Access polygon Points from self.inlet.polygon 305 306 # The Following computes the flow crossing over 3 sides of the exchange polygon for the structure 307 # Clearly the flow in the culvert can not be more than that flowing toward it through the exhange polygon 308 309 #q1 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][1:3]) # First Side Segment 310 #q2 = domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'][2:]) # Second Face Segment 311 #q3 =domain.get_flow_through_cross_section(self.culvert_polygons['exchange_polygon0'].take([3,0], axis=0)) # Third Side Segment 312 # q4 = domain.get_flow_through_cross_section([self.culvert_polygons['exchange_polygon0'][1:4]][0]) 313 #q4=max(q1,0.0)+max(q2,0.0)+max(q3,0.0) 314 # To use only the Flow crossing the 3 sides of the Exchange Polygon use the following Line Only 315 #max_Q=max(q1,q2,q3,q4) 316 # Try Simple Smoothing using Average of 2 approaches 317 #max_Q=(max(q1,q2,q3,q4)+max_Q)/2.0 318 # Calculate the minimum in absolute terms of 319 # the requsted flow and the possible flow 320 Q_reduced = sign(Q)*min(abs(Q), abs(max_Q)) 321 if self.verbose is True: 322 msg = 'Initial Q Reduced = %.2f m3/s. ' % (Q_reduced) 323 log.critical(msg) 324 325 if self.log_filename is not None: 326 log_to_file(self.log_filename, msg) 327 # Now Keep Rolling Average of Computed Discharge to Reduce / Remove Oscillations 328 # can use delta_t if we want to averageover a time frame for example 329 # N = 5.0/delta_t Will provide the average over 5 seconds 330 331 self.i=(self.i+1)%self.N 332 self.Q_list[self.i]=Q_reduced 333 Q_reduced = sum(self.Q_list)/len(self.Q_list) 334 335 if self.verbose is True: 336 msg = 'Final Q Reduced = %.2f m3/s. ' % (Q_reduced) 337 log.critical(msg) 338 339 if self.log_filename is not None: 340 log_to_file(self.log_filename, msg) 341 342 343 if abs(Q_reduced) < abs(Q): 344 msg = '%.2fs: Requested flow is ' % time 345 msg += 'greater than what is supported by the smallest ' 346 msg += 'depth at inlet exchange area:\n ' 347 msg += 'inlet exchange area: %.2f '% (I.exchange_area) 348 msg += 'velocity at inlet :%.2f '% (I.velocity) 349 msg += 'Vel* Exch Area = : %.2f '% (I.velocity*avg_depth*self.width) 350 msg += 'h_min*inlet_area/delta_t = %.2f*%.2f/%.2f '\ 351 % (avg_depth, I.exchange_area, delta_t) 352 msg += ' = %.2f m^3/s\n ' % Q_reduced 353 msg += 'Q will be reduced from %.2f m^3/s to %.2f m^3/s.' % (Q, Q_reduced) 354 msg += 'Note calculate max_Q from V %.2f m^3/s ' % (max_Q) 355 if self.verbose is True: 356 log.critical(msg) 357 358 if self.log_filename is not None: 359 log_to_file(self.log_filename, msg) 360 361 return Q_reduced 233 362 234 363 364 def compute_rates(self, delta_t):365 """Compute new rates for inlet and outlet366 """367 # Time stuff368 time = self.domain.get_time()369 self.last_update = time370 371 if hasattr(self, 'log_filename'):372 log_filename = self.log_filename373 374 # Compute stage, energy and velocity at the375 # enquiry points at each end of the culvert376 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_velocities387 else:388 u = 0.0389 v = 0.0390 391 uv2 = u**2 + v**2392 393 if self.use_velocity_head is True:394 velocity_head = 0.5*uv2/g395 else:396 velocity_head = 0.0397 398 inlet.set_total_energy(velocity_head + stage)399 inlet.set_specific_energy(velocity_head + depth)400 401 # We now need to deal with each opening individually402 # Determine flow direction based on total energy difference403 delta_total_energy = inlets[0].get_total_energy() - inlets[1].get_total_energy()404 405 if delta_total_energy >= 0:406 inlet = inlets[0]407 outlet = inlets[1]408 else:409 inlet = inlets[1]410 outlet = inlets[2]411 412 #msg = 'Total energy difference is negative'413 #assert delta_total_energy >= 0.0, msg414 415 # Recompute slope and issue warning if flow is uphill416 # These values do not enter the computation417 delta_z = inlet.get_average_elevation() - outlet.get_average_elevation()418 culvert_slope = (delta_z/self.culvert_length)419 if culvert_slope < 0.0:420 # Adverse gradient - flow is running uphill421 # Flow will be purely controlled by uphill outlet face422 if self.verbose is True:423 log.critical('%.2fs - WARNING: Flow is running uphill.' % time)424 425 if self.log_filename is not None:426 s = 'Time=%.2f, inlet stage = %.2f, outlet stage = %.2f'\427 %(time, self.inlet.stage, self.outlet.stage)428 log_to_file(self.log_filename, s)429 s = 'Delta total energy = %.3f' %(delta_total_energy)430 log_to_file(log_filename, s)431 432 # Determine controlling energy (driving head) for culvert433 if inlet.get_specific_energy() > delta_total_energy:434 # Outlet control435 driving_head = delta_total_energy436 else:437 # Inlet control438 driving_head = inlet.get_specific_energy()439 440 inlet_depth = inlet.get_average_stage() - inlet.get_average_elevation()441 outlet_depth = outlet.get_average_stage() - outlet.get_average_elevation()442 443 if inlet_depth <= self.trigger_depth:444 Q = 0.0445 else:446 # Calculate discharge for one barrel and447 # set inlet.rate and outlet.rate448 449 if self.culvert_description_filename is not None:450 try:451 Q = interpolate_linearly(driving_head,452 self.rating_curve[:,0],453 self.rating_curve[:,1])454 except Below_interval, e:455 Q = self.rating_curve[0,1]456 msg = '%.2fs: ' % time457 msg += 'Delta head smaller than rating curve minimum: '458 msg += str(e)459 msg += '\n '460 msg += 'I will use minimum discharge %.2f m^3/s ' % Q461 msg += 'for culvert "%s"' % self.label462 463 if hasattr(self, 'log_filename'):464 log_to_file(self.log_filename, msg)465 except Above_interval, e:466 Q = self.rating_curve[-1,1]467 msg = '%.2fs: ' % time468 msg += 'Delta head greater than rating curve maximum: '469 msg += str(e)470 msg += '\n '471 msg += 'I will use maximum discharge %.2f m^3/s ' % Q472 msg += 'for culvert "%s"' % self.label473 474 if self.log_filename is not None:475 log_to_file(self.log_filename, msg)476 else:477 # User culvert routine478 Q, barrel_velocity, culvert_outlet_depth =\479 self.culvert_routine(inlet_depth,480 outlet_depth,481 inlet.get_average_velocity(),482 outlet.get_average_velocity(),483 inlet.get_specific_energy(),484 delta_total_energy,485 g,486 culvert_length=self.length,487 culvert_width=self.width,488 culvert_height=self.height,489 culvert_type=self.culvert_type,490 manning=self.manning,491 sum_loss=self.sum_loss,492 log_filename=self.log_filename)493 494 # Adjust discharge for multiple barrels495 Q *= self.number_of_barrels496 497 # Adjust discharge for available water at the inlet498 Q = self.adjust_flow_for_available_water_at_inlet(Q, delta_t)499 500 self.inlet.rate = -Q501 self.outlet.rate = Q502 503 # Momentum jet stuff504 if self.use_momentum_jet is True:505 506 # Compute barrel momentum507 barrel_momentum = barrel_velocity*culvert_outlet_depth508 509 if self.log_filename is not None:510 s = 'Barrel velocity = %f' %barrel_velocity511 log_to_file(self.log_filename, s)512 513 # Compute momentum vector at outlet514 outlet_mom_x, outlet_mom_y = self.vector * barrel_momentum515 516 if self.log_filename is not None:517 s = 'Directional momentum = (%f, %f)' %(outlet_mom_x, outlet_mom_y)518 log_to_file(self.log_filename, s)519 520 # Update momentum521 if delta_t > 0.0:522 xmomentum_rate = outlet_mom_x - outlet.momentum[0].value523 xmomentum_rate /= delta_t524 525 ymomentum_rate = outlet_mom_y - outlet.momentum[1].value526 ymomentum_rate /= delta_t527 528 if self.log_filename is not None:529 s = 'X Y MOM_RATE = (%f, %f) ' %(xmomentum_rate, ymomentum_rate)530 log_to_file(self.log_filename, s)531 else:532 xmomentum_rate = ymomentum_rate = 0.0533 534 535 # Set momentum rates for outlet jet536 outlet.momentum[0].rate = xmomentum_rate537 outlet.momentum[1].rate = ymomentum_rate538 539 # Remember this value for next step (IMPORTANT)540 outlet.momentum[0].value = outlet_mom_x541 outlet.momentum[1].value = outlet_mom_y542 543 if int(domain.time*100) % 100 == 0:544 545 if self.log_filename is not None:546 s = 'T=%.5f, Culvert Discharge = %.3f f'\547 %(time, Q)548 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\549 %(culvert_outlet_depth, outlet_mom_x,outlet_mom_y)550 s += ' Momentum rate: (%.4f, %.4f)'\551 %(xmomentum_rate, ymomentum_rate)552 s+='Outlet Vel= %.3f'\553 %(barrel_velocity)554 log_to_file(self.log_filename, s)555 556 557 # Execute momentum terms558 # This is where Inflow objects are evaluated and update the domain559 self.outlet.momentum[0](domain)560 self.outlet.momentum[1](domain)561 562 563 564 # Log timeseries to file565 try:566 fid = open(self.timeseries_filename, 'a')567 except:568 pass569 else:570 fid.write('%.2f, %.2f\n' %(time, Q))571 fid.close()572 573 # Store value of time574 self.last_time = time575 576 577 235 # FIXME(Ole): Write in C and reuse this function by similar code 578 236 # in interpolate.py -
trunk/anuga_core/source/anuga/structures/inlet.py
r7974 r7975 1 1 2 from anuga.geometry.polygon import inside_polygon, is_inside_polygon 2 from anuga.config import velocity_protection 3 from anuga.config import velocity_protection, g 3 4 import math 4 5 … … 117 118 def get_total_water_volume(self): 118 119 119 return num.sum(self.get_heights *self.get_areas())120 return num.sum(self.get_heights()*self.get_areas()) 120 121 121 122 … … 134 135 135 136 136 def get_average_ velocity(self):137 def get_average_speed(self): 137 138 138 139 u, v = self.get_velocities() … … 141 142 average_v = num.sum(v)/self.triangle_indices.size 142 143 143 return math.sqrt(average_u**2 + average_v**2) 144 return math.sqrt(average_u**2 + average_v**2) 144 145 145 146 146 def set_total_energy(self, total_energy): 147 148 def get_average_velocity_head(self): 149 150 return 0.5*self.get_average_speed()**2/g 151 152 153 def get_average_total_energy(self): 147 154 148 self.total_energy = total_energy155 return self.get_average_velocity_head() + self.get_average_stage() 149 156 157 158 def get_average_specific_energy(self): 150 159 151 def get_total_energy(self): 152 153 return self.total_energy 154 155 156 def set_specific_energy(self, specific_energy): 157 158 self.specific_energy = specific_energy 160 return self.get_average_velocity_head() + self.get_average_height() 159 161 160 161 def get_specific_energy(self): 162 163 return self.specific_energy 164 165 162 163 164 def set_heights(self,height): 165 166 self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, self.get_elevations() + height) 167 168 169 def set_stages(self,stage): 170 171 self.domain.quantities['stage'].centroid_values.put(self.triangle_indices, stage) 172 173 def set_xmoms(self,xmom): 174 175 self.xmoms=self.domain.quantities['xmomentum'].centroid_values.put(self.triangle_indices, xmom) 176 177 178 def set_ymoms(self,ymom): 179 180 self.xmoms=self.domain.quantities['ymomentum'].centroid_values.put(self.triangle_indices, ymom) 181 182 183 def set_elevations(self,elevation): 184 185 self.domain.quantities['elevation'].centroid_values.put(self.triangle_indices, elevation) -
trunk/anuga_core/source/anuga/structures/testing_culvert.py
r7968 r7975 42 42 domain.set_name('Test_culvert') # Output name 43 43 domain.set_default_order(2) 44 #domain.set_beta(1.5) 44 45 45 46 … … 110 111 111 112 ## Inflow based on Flow Depth and Approaching Momentum 112 Bi = anuga.Dirichlet_boundary([ 1.0, 0.0, 0.0])113 Bi = anuga.Dirichlet_boundary([2.0, 0.0, 0.0]) 113 114 Br = anuga.Reflective_boundary(domain) # Solid reflective wall 114 115 #Bo = anuga.Dirichlet_boundary([-5, 0, 0]) # Outflow … … 131 132 #min_delta_w = sys.maxint 132 133 #max_delta_w = -min_delta_w 133 for t in domain.evolve(yieldstep = 1.0, finaltime = 100):134 for t in domain.evolve(yieldstep = 1.0, finaltime = 200): 134 135 domain.write_time() 136 137 if domain.get_time() > 150.5 and domain.get_time() < 151.5 : 138 Bi = anuga.Dirichlet_boundary([0.0, 0.0, 0.0]) 139 domain.set_boundary({'left': Bi, 'right': Br, 'top': Br, 'bottom': Br}) 140 135 141 #delta_w = culvert.inlet.stage - culvert.outlet.stage 136 142
Note: See TracChangeset
for help on using the changeset viewer.
