 Timestamp:
 May 16, 2008, 5:12:37 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

anuga_work/development/Rudy_vandrie_2007/culvert_development.py
r5314 r5340 286 286 openings = self.openings 287 287 for i, opening in enumerate(openings): 288 stage = domain.quantities['stage'].get_values(location='centroids', 288 # Quantities corresponding to fluid exchange field for this opening 289 stage_o = domain.quantities['stage'].get_values(location='centroids', 289 290 indices=opening.indices) 290 elevation = domain.quantities['elevation'].get_values(location='centroids',291 elevation_o = domain.quantities['elevation'].get_values(location='centroids', 291 292 indices=opening.indices) 293 xmomentum_o = domain.quantities['xmomentum'].get_values(location='centroids', 294 indices=opening.indices) 295 ymomentum_o = domain.quantities['xmomentum'].get_values(location='centroids', 296 indices=opening.indices) 297 298 # Compute mean velocity in the exchange area in front of the culvert (taking zero depths into account) 299 depth_o = stage_o  elevation_o 300 301 ux_o = xmomentum_o/(depth_o+epsilon) 302 uy_o = ymomentum_o/(depth_o+epsilon) 303 v_o = mean(sqrt(ux_o**2+uy_o**2)) 304 d_o = mean(depth_o) 305 w_o = mean(stage_o) 306 z_o = mean(elevation_o) 307 self.mean_xmomentum_o = mean(xmomentum_o) 308 self.mean_ymomentum_o = mean(ymomentum_o) 292 309 293 310 # Indices corresponding to energy enquiry field for this opening … … 317 334 w = mean(stage) 318 335 z = mean(elevation) 319 336 self.mean_xmomentum = mean(xmomentum) 337 self.mean_ymomentum = mean(ymomentum) 338 320 339 # Ratio of depth to culvert height 321 340 ratio = d/(2*self.height) … … 436 455 delta_t = time  self.last_time 437 456 if delta_t > 0.0: 438 xmomentum_rate = outlet_mom_x  self.xmom_forcing1.value 439 xmomentum_rate /= delta_t 457 xmomentum_rate = outlet_mom_x  self.mean_xmomentum 458 if xmomentum_rate > 0: 459 xmomentum_rate /= delta_t 460 else: 461 xmomentum_rate = 0.0 440 462 441 ymomentum_rate = outlet_mom_y  self.ymom_forcing1.value 442 ymomentum_rate /= delta_t 463 ymomentum_rate = outlet_mom_y  self.mean_ymomentum 464 if ymomentum_rate > 0: 465 ymomentum_rate /= delta_t 466 else: 467 ymomentum_rate = 0.0 443 468 else: 444 469 xmomentum_rate = ymomentum_rate = 0.0 … … 453 478 454 479 455 if int( domain.time*100) % 100 == 0:456 s = 'T=%. 2f, Culvert Discharge = %.3f Culv. Vel. %0.3f'\480 if int(time*100) % 100 == 0: 481 s = 'T=%.3f, Culvert Discharge = %.3f Culv. Vel. %0.3f'\ 457 482 %(time, flow_rate_control, outlet_vel) 458 s += ' Depth= %0.3f Momentum = (%0.3f, %0.3f)'\ 459 %(outlet_dep, outlet_mom_x,outlet_mom_y) 460 s += ' Momentum rate: (%.4f, %.4f)'\ 483 s += ' Depth= %0.3f\n Outlet Momentum = (%0.3f, %0.3f)\n'\ 484 %(outlet_dep, outlet_mom_x, outlet_mom_y) 485 s += ' Avg Momentum at opening = (%0.3f, %0.3f)\n'\ 486 %(self.mean_xmomentum_o, self.mean_ymomentum_o) 487 s += ' Avg Momentum in enquiry = (%0.3f, %0.3f)\n'\ 488 %(self.mean_xmomentum, self.mean_ymomentum) 489 s += ' Momentum rate: (%.4f, %.4f)'\ 461 490 %(xmomentum_rate, ymomentum_rate) 462 491 fid.write(s) … … 488 517 # Print out flows every 1 seconds 489 518 if int(time*100) % 100 == 0: 490 s = 'Time %. 2f\n' %time519 s = 'Time %.3f\n' %time 491 520 s += ' Opening[0]: d=%.3f, A=%.3f, E=%.3f, r=%.3f\n'\ 492 521 %(openings[0].depth, … … 506 535 507 536 print s 537 508 538 fid.write(s) 509 539
Note: See TracChangeset
for help on using the changeset viewer.