Ignore:
Timestamp:
Oct 5, 2006, 5:50:11 PM (17 years ago)
Author:
steve
Message:
 
File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/shallow_water_domain.py

    r3689 r3703  
    144144        self.beta_h      = beta_h
    145145
    146 
     146        self.flux_function = flux_function_central
     147        #self.flux_function = flux_function_kinetic
     148       
    147149        self.forcing_terms.append(manning_friction)
    148150        self.forcing_terms.append(gravity)
     
    495497####################################
    496498# Flux computation
    497 def flux_function(normal, ql, qr, zl, zr):
     499def flux_function_central(normal, ql, qr, zl, zr):
    498500    """Compute fluxes between volumes for the shallow water wave equation
    499501    cast in terms of w = h+z using the 'central scheme' as described in
     
    577579
    578580    return edgeflux, max_speed
     581
     582def 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
     596def 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
    579686
    580687
     
    648755
    649756            #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)
    651758            flux -= edgeflux * domain.edgelengths[k,i]
    652759
     
    712819
    713820    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,
    716823                                     domain.neighbours,
    717824                                     domain.neighbour_edges,
Note: See TracChangeset for help on using the changeset viewer.