- Timestamp:
- Oct 5, 2006, 5:50:11 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/shallow_water_domain.py
r3689 r3703 144 144 self.beta_h = beta_h 145 145 146 146 self.flux_function = flux_function_central 147 #self.flux_function = flux_function_kinetic 148 147 149 self.forcing_terms.append(manning_friction) 148 150 self.forcing_terms.append(gravity) … … 495 497 #################################### 496 498 # Flux computation 497 def flux_function (normal, ql, qr, zl, zr):499 def flux_function_central(normal, ql, qr, zl, zr): 498 500 """Compute fluxes between volumes for the shallow water wave equation 499 501 cast in terms of w = h+z using the 'central scheme' as described in … … 577 579 578 580 return edgeflux, max_speed 581 582 def erfcc(x): 583 584 from math import fabs, exp 585 586 z=fabs(x) 587 t=1.0/(1.0+0.5*z) 588 result=t*exp(-z*z-1.26551223+t*(1.00002368+t*(.37409196+ 589 t*(.09678418+t*(-.18628806+t*(.27886807+t*(-1.13520398+ 590 t*(1.48851587+t*(-.82215223+t*.17087277))))))))) 591 if x < 0.0: 592 result = 2.0-result 593 594 return result 595 596 def flux_function_kinetic(normal, ql, qr, zl, zr): 597 """Compute fluxes between volumes for the shallow water wave equation 598 cast in terms of w = h+z using the 'central scheme' as described in 599 600 Zhang et. al., Advances in Water Resources, 26(6), 2003, 635-647. 601 602 603 Conserved quantities w, uh, vh are stored as elements 0, 1 and 2 604 in the numerical vectors ql an qr. 605 606 Bed elevations zl and zr. 607 """ 608 609 from anuga.config import g, epsilon 610 from math import sqrt 611 from Numeric import array 612 613 #Align momentums with x-axis 614 q_left = rotate(ql, normal, direction = 1) 615 q_right = rotate(qr, normal, direction = 1) 616 617 z = (zl+zr)/2 #Take average of field values 618 619 w_left = q_left[0] #w=h+z 620 h_left = w_left-z 621 uh_left = q_left[1] 622 623 if h_left < epsilon: 624 u_left = 0.0 #Could have been negative 625 h_left = 0.0 626 else: 627 u_left = uh_left/h_left 628 629 630 w_right = q_right[0] #w=h+z 631 h_right = w_right-z 632 uh_right = q_right[1] 633 634 635 if h_right < epsilon: 636 u_right = 0.0 #Could have been negative 637 h_right = 0.0 638 else: 639 u_right = uh_right/h_right 640 641 vh_left = q_left[2] 642 vh_right = q_right[2] 643 644 soundspeed_left = sqrt(g*h_left) 645 soundspeed_right = sqrt(g*h_right) 646 647 #Maximal wave speed 648 s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0) 649 650 #Minimal wave speed 651 s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0) 652 653 #Flux computation 654 655 F_left = 0.0 656 F_right = 0.0 657 from math import sqrt, pi, exp 658 if h_left > 0.0: 659 F_left = u_left/sqrt(g*h_left) 660 if h_right > 0.0: 661 F_right = u_right/sqrt(g*h_right) 662 663 edgeflux = array([0.0, 0.0, 0.0]) 664 665 edgeflux[0] = h_left*u_left/2.0*erfcc(-F_left) + \ 666 h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ 667 h_right*u_right/2.0*erfcc(F_right) - \ 668 h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) 669 670 edgeflux[1] = (h_left*u_left**2 + g/2.0*h_left**2)/2.0*erfcc(-F_left) + \ 671 u_left*h_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ 672 (h_right*u_right**2 + g/2.0*h_right**2)/2.0*erfcc(F_right) - \ 673 u_right*h_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) 674 675 edgeflux[2] = vh_left*u_left/2.0*erfcc(-F_left) + \ 676 vh_left*sqrt(g*h_left)/2.0/sqrt(pi)*exp(-(F_left**2)) + \ 677 vh_right*u_right/2.0*erfcc(F_right) - \ 678 vh_right*sqrt(g*h_right)/2.0/sqrt(pi)*exp(-(F_right**2)) 679 680 681 edgeflux = rotate(edgeflux, normal, direction=-1) 682 max_speed = max(abs(s_max), abs(s_min)) 683 684 return edgeflux, max_speed 685 579 686 580 687 … … 648 755 649 756 #Flux computation using provided function 650 edgeflux, max_speed = flux_function(normal, ql, qr, zl, zr)757 edgeflux, max_speed = domain.flux_function(normal, ql, qr, zl, zr) 651 758 flux -= edgeflux * domain.edgelengths[k,i] 652 759 … … 712 819 713 820 timestep = float(sys.maxint) 714 from shallow_water_ext import compute_fluxes 715 domain.timestep = compute_fluxes (timestep, domain.epsilon, domain.g,821 from shallow_water_ext import compute_fluxes_ext_central as compute_fluxes_ext 822 domain.timestep = compute_fluxes_ext(timestep, domain.epsilon, domain.g, 716 823 domain.neighbours, 717 824 domain.neighbour_edges,
Note: See TracChangeset
for help on using the changeset viewer.