Changeset 232 for inundation/ga/storm_surge/pyvolution/shallow_water.py

Timestamp:

Aug 27, 2004, 1:35:01 PM (20 years ago)

Author:

ole

Message:

Fixed flux_function to match c-version
More testing

File:

• inundation/ga/storm_surge/pyvolution/shallow_water.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/shallow_water.py

@@ -98 +98 @@
     """
 
-    from config import g
+    from config import g, epsilon
     from math import sqrt
     from Numeric import array
@@ -114 +114 @@
     uh_left = q_left[1]
 
-    try:
+    if h_left < epsilon:
+        u_left = 0.0  #Could have been negative
+        h_left = 0.0
+    else:    
         u_left  = uh_left/h_left
-    except:
-        u_left = 0.0
-        h_left = 0.0
 
 
@@ -126 +126 @@
 
 
-    try:
+    if h_right < epsilon:
+        u_right = 0.0  #Could have been negative
+        h_right = 0.0
+    else:    
         u_right  = uh_right/h_right
-    except:
-        u_right = 0.0
-        h_right = 0.0
-
-        
+
     vh_left  = q_left[2]
     vh_right = q_right[2]        
@@ -146 +145 @@
 
     #Flux computation
-    flux_left  = array([uh_left,
-                        u_left**2*h_left + 0.5*g*h_left**2,
-                        vh_left*u_left])
-    flux_right = array([uh_right,
-                        u_right**2*h_right + 0.5*g*h_right**2,
-                        vh_right*u_right])    
+    flux_left  = array([u_left*h_left,
+                        u_left*uh_left + 0.5*g*h_left**2,
+                        u_left*vh_left])
+    flux_right = array([u_right*h_right,
+                        u_right*uh_right + 0.5*g*h_right**2,
+                        u_right*vh_right])    
 
     denom = s_max-s_min
     if denom == 0.0:
-        flux = array([0.0, 0.0, 0.0])
+        edgeflux = array([0.0, 0.0, 0.0])
         max_speed = 0.0
     else:    
-        flux = (s_max*flux_left - s_min*flux_right)/denom
-        flux += s_max*s_min*(q_right-q_left)/denom
-
-        flux = rotate(flux, normal, direction=-1)
+        edgeflux = (s_max*flux_left - s_min*flux_right)/denom
+        edgeflux += s_max*s_min*(q_right-q_left)/denom
+
+        edgeflux = rotate(edgeflux, normal, direction=-1)
         max_speed = max(abs(s_max), abs(s_min))
 
-    return flux, max_speed        
+    return edgeflux, max_speed