Changeset 5563 for anuga_work/development/anuga_1d/test_shallow_water.py

Timestamp:

Jul 23, 2008, 4:38:10 PM (15 years ago)

Author:

steve

Message:

Working version of comp_flux_ext

File:

• anuga_work/development/anuga_1d/test_shallow_water.py (modified) (3 diffs)

anuga_work/development/anuga_1d/test_shallow_water.py

@@ -41 +41 @@
         assert allclose( domain.quantities['stage'].explicit_update, stage_ud )
         assert allclose( domain.quantities['xmomentum'].explicit_update, xmom_ud )
+
+
+    def test_local_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        assert allclose(array([2.0, 8.9],Float), edgeflux)
+        assert allclose(5.1304951685, maxspeed)
+
+        normal = -1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        assert allclose(array([-2.0, -8.9],Float), edgeflux)
+        assert allclose(5.1304951685, maxspeed)
+
+    def test_domain_flux_function(self):
+        normal = 1.0
+        ql = array([1.0, 2.0],Float)
+        qr = array([1.0, 2.0],Float)
+        zl = 0.0
+        zr = 0.0
+
+        edgeflux, maxspeed = flux_function(normal, ql,qr,zl,zr)
+
+        #print edgeflux
+
+        from shallow_water_domain import flux_function as domain_flux_function
+
+        domainedgeflux, domainmaxspeed = domain_flux_function(normal, ql,qr,zl,zr)
+
+        #print domainedgeflux
+
+        assert allclose(domainedgeflux, edgeflux)
+        assert allclose(domainmaxspeed, maxspeed)
+
 
 
@@ -128 +173 @@
                 #m = domain.neighbour_edges[k,i]
                 m = domain.neighbour_vertices[k,i]
+                #print i, ' ' , m
                 #qr = [stage[n, m], xmom[n, m], ymom[n, m]]
                 qr[0] = stage[n, m]
@@ -169 +215 @@
 
 
+def flux_function(normal, ql, qr, zl, zr):
+    """Compute fluxes between volumes for the shallow water wave equation
+    cast in terms of w = h+z using the 'central scheme' as described in
+
+    Kurganov, Noelle, Petrova. 'Semidiscrete Central-Upwind Schemes For
+    Hyperbolic Conservation Laws and Hamilton-Jacobi Equations'.
+    Siam J. Sci. Comput. Vol. 23, No. 3, pp. 707-740.
+
+    The implemented formula is given in equation (3.15) on page 714
+
+    Conserved quantities w, uh, are stored as elements 0 and 1
+    in the numerical vectors ql an qr.
+
+    Bed elevations zl and zr.
+    """
+
+    from config import g, epsilon
+    from math import sqrt
+    from Numeric import array
+
+    #print 'ql',ql
+
+    #Align momentums with x-axis
+    #q_left  = rotate(ql, normal, direction = 1)
+    #q_right = rotate(qr, normal, direction = 1)
+    q_left = ql
+    q_left[1] = q_left[1]*normal
+    q_right = qr
+    q_right[1] = q_right[1]*normal
+
+    #z = (zl+zr)/2 #Take average of field values
+    z = 0.5*(zl+zr) #Take average of field values
+
+    w_left  = q_left[0]   #w=h+z
+    h_left  = w_left-z
+    uh_left = q_left[1]
+
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:
+        u_left  = uh_left/h_left
+
+
+    w_right  = q_right[0]  #w=h+z
+    h_right  = w_right-z
+    uh_right = q_right[1]
+
+
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:
+        u_right  = uh_right/h_right
+
+    #vh_left  = q_left[2]
+    #vh_right = q_right[2]
+
+    #print h_right
+    #print u_right
+    #print h_left
+    #print u_right
+
+    soundspeed_left  = sqrt(g*h_left)
+    soundspeed_right = sqrt(g*h_right)
+
+    #Maximal wave speed
+    s_max = max(u_left+soundspeed_left, u_right+soundspeed_right, 0)
+    
+    #Minimal wave speed
+    s_min = min(u_left-soundspeed_left, u_right-soundspeed_right, 0)
+    
+    #Flux computation
+
+    #flux_left  = array([u_left*h_left,
+    #                    u_left*uh_left + 0.5*g*h_left**2])
+    #flux_right = array([u_right*h_right,
+    #                    u_right*uh_right + 0.5*g*h_right**2])
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left*h_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right*h_right])
+
+    denom = s_max-s_min
+    if denom == 0.0:
+        edgeflux = array([0.0, 0.0])
+        max_speed = 0.0
+    else:
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+        
+        edgeflux[1] = edgeflux[1]*normal
+
+        max_speed = max(abs(s_max), abs(s_min))
+
+    return edgeflux, max_speed
+
 
 #-------------------------------------------------------------