Changeset 5536 for anuga_work/development/anuga_1d/shallow_water_domain.py
- Timestamp:
- Jul 18, 2008, 4:37:03 PM (16 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_work/development/anuga_1d/shallow_water_domain.py
r5535 r5536 45 45 #from domain import * 46 46 #from domain_order2 import * 47 from domain _t2import *47 from domain import * 48 48 Generic_Domain = Domain #Rename 49 49 … … 388 388 389 389 390 def flux_function_split(normal, ql, qr, zl, zr): 391 from config import g, epsilon 392 from math import sqrt 393 from Numeric import array 394 395 #print 'ql',ql 396 397 #Align momentums with x-axis 398 #q_left = rotate(ql, normal, direction = 1) 399 #q_right = rotate(qr, normal, direction = 1) 400 q_left = ql 401 q_left[1] = q_left[1]*normal 402 q_right = qr 403 q_right[1] = q_right[1]*normal 404 405 #z = (zl+zr)/2 #Take average of field values 406 z = 0.5*(zl+zr) #Take average of field values 407 408 w_left = q_left[0] #w=h+z 409 h_left = w_left-z 410 uh_left = q_left[1] 411 412 if h_left < epsilon: 413 u_left = 0.0 #Could have been negative 414 h_left = 0.0 415 else: 416 u_left = uh_left/h_left 417 418 419 w_right = q_right[0] #w=h+z 420 h_right = w_right-z 421 uh_right = q_right[1] 422 423 424 if h_right < epsilon: 425 u_right = 0.0 #Could have been negative 426 h_right = 0.0 427 else: 428 u_right = uh_right/h_right 429 430 #vh_left = q_left[2] 431 #vh_right = q_right[2] 432 433 #soundspeed_left = sqrt(g*h_left) 434 #soundspeed_right = sqrt(g*h_right) 435 436 #Maximal wave speed 437 #s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0) 438 s_max = max(u_left, u_right, 0) 439 440 #Minimal wave speed 441 #s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0) 442 s_min = min(u_left, u_right, 0) 443 444 #Flux computation 445 446 #flux_left = array([u_left*h_left, 447 # u_left*uh_left + 0.5*g*h_left*h_left]) 448 #flux_right = array([u_right*h_right, 449 # u_right*uh_right + 0.5*g*h_right*h_right]) 450 flux_left = array([u_left*h_left, 451 u_left*uh_left])# + 0.5*g*h_left*h_left]) 452 flux_right = array([u_right*h_right, 453 u_right*uh_right])# + 0.5*g*h_right*h_right]) 454 455 denom = s_max-s_min 456 if denom == 0.0: 457 edgeflux = array([0.0, 0.0]) 458 max_speed = 0.0 459 else: 460 edgeflux = (s_max*flux_left - s_min*flux_right)/denom 461 edgeflux += s_max*s_min*(q_right-q_left)/denom 462 463 edgeflux[1] = edgeflux[1]*normal 464 465 max_speed = max(abs(s_max), abs(s_min)) 466 467 return edgeflux, max_speed 390 468 391 469 392 def compute_timestep(domain): … … 518 441 normal = domain.normals[k, i] 519 442 520 if domain.split == False: 521 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr) 522 elif domain.split == True: 523 edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr) 443 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr) 444 524 445 #Update optimal_timestep 525 446 try: … … 612 533 #qr = [stage[n, m], xmom[n, m], ymom[n, m]] 613 534 qr[0] = stage[n, m] 614 qr[1] = 535 qr[1] = xmom[n, m] 615 536 zr = bed[n, m] 616 537 … … 625 546 626 547 627 if domain.split == False: 628 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr) 629 elif domain.split == True: 630 edgeflux, max_speed = flux_function_split(normal, ql, qr, zl, zr) 548 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr) 549 631 550 #print 'edgeflux', edgeflux 632 551 … … 646 565 647 566 Stage.explicit_update[k] = flux[0] 648 Xmom.explicit_update[k] = flux[1]567 Xmom.explicit_update[k] = flux[1] 649 568 #Ymom.explicit_update[k] = flux[2] 650 569 #print "flux cell",k,flux[0] … … 1054 973 1055 974 #Update momentum (explicit update is reset to source values) 1056 if domain.split == False: 1057 xmom[k] += -g*bx*avg_h 1058 #xmom[k] = -g*bx*avg_h 1059 #stage[k] = 0.0 1060 elif domain.split == True: 1061 xmom[k] += -g*wx*avg_h 1062 #xmom[k] = -g*wx*avg_h 1063 #ymom[k] += -g*zy*avg_h 975 xmom[k] += -g*bx*avg_h 976 #xmom[k] = -g*bx*avg_h 977 #stage[k] = 0.0 978 1064 979 1065 980 def manning_friction(domain):
Note: See TracChangeset
for help on using the changeset viewer.