Changeset 8374 for trunk/anuga_work/shallow_water_balanced_steve

Timestamp:

Mar 30, 2012, 5:14:21 PM (13 years ago)

Author:

steve

Message:

Still working on channel flow

Location:

trunk/anuga_work/shallow_water_balanced_steve

Files:

• swb_domain.py (modified) (6 diffs)
• swb_domain_ext.c (modified) (1 diff)

trunk/anuga_work/shallow_water_balanced_steve/swb_domain.py

@@ -130 +130 @@
 
 
-    def distribute_to_vertices_and_edges(self):
+    def distribute_to_vertices_and_edges_w_v(self):
         """Distribution from centroids to edges specific to the SWW eqn.
 
@@ -203 +203 @@
                 raise Exception('Unknown order: %s' % str(self._order_))
 
-        for name in [ 'height' ]:
-            Q = self.quantities[name]
-            if self._order_ == 1:
-                Q.extrapolate_first_order()
-            elif self._order_ == 2:
-                #Q.extrapolate_second_order_and_limit_by_edge()
-                Q.extrapolate_second_order_and_limit_by_vertex()
-            else:
-                raise Exception('Unknown order: %s' % str(self._order_))
 
 
@@ -227 +218 @@
         #protect(self)
 
-        z_V[:]  = w_V - h_V
+        h_V[:]  = w_V - z_V
         uh_V[:] = u_V * h_V
         vh_V[:] = v_V * h_V
@@ -246 +237 @@
 
 
-
-    def distribute_to_vertices_and_edges_h(self):
+    def distribute_to_vertices_and_edges_w_h(self):
         """Distribution from centroids to edges specific to the SWW eqn.
 
@@ -256 +246 @@
         In addition, all conserved quantities get distributed as per either a
         constant (order==1) or a piecewise linear function (order==2).
-        
+
 
         Precondition:
         All conserved quantities defined at centroids and bed elevation defined at
         edges.
-        
+
         Postcondition
         Evolved quantities defined at vertices and edges
         """
 
+        from anuga.shallow_water.shallow_water_domain import \
+                  protect_against_infinitesimal_and_negative_heights as protect
 
         #Shortcuts
@@ -276 +268 @@
         V  = self.quantities['yvelocity']
 
+        #Arrays
+        w_C   = W.centroid_values
+        uh_C  = UH.centroid_values
+        vh_C  = VH.centroid_values
+        z_C   = Z.centroid_values
+        h_C   = H.centroid_values
+        u_C   = U.centroid_values
+        v_C   = V.centroid_values
+
+        #num_min = num.min(w_C-z_C)
+        #if num_min < -1.0e-5:
+        #    print '**** num.min(w_C-z_C)', num_min
+
+
+        w_C[:] = num.maximum(w_C, z_C)
+        h_C[:] = w_C - z_C
+
+
+        assert num.min(h_C) >= 0.0
+
+        num.putmask(uh_C, h_C < 1.0e-15, 0.0)
+        num.putmask(vh_C, h_C < 1.0e-15, 0.0)
+        #num.putmask(h_C, h_C < 1.0e-15, 1.0e-16)
+
+        # Noelle has an alternative method for calculating velcities
+        # Check it out in the GPU shallow water paper
+        H0 = 1.0e-16
+        u_C[:]  = uh_C/(h_C + H0/h_C)
+        v_C[:]  = vh_C/(h_C + H0/h_C)
+
+        #num.putmask(h_C, h_C < 1.0e-15, 0.0)
+
+        for name in [ 'stage', 'xvelocity', 'yvelocity' ]:
+            Q = self.quantities[name]
+            if self._order_ == 1:
+                Q.extrapolate_first_order()
+            elif self._order_ == 2:
+                Q.extrapolate_second_order_and_limit_by_edge()
+                #Q.extrapolate_second_order_and_limit_by_vertex()
+            else:
+                raise Exception('Unknown order: %s' % str(self._order_))
+
+        for name in [ 'height' ]:
+            Q = self.quantities[name]
+            if self._order_ == 1:
+                Q.extrapolate_first_order()
+            elif self._order_ == 2:
+                #Q.extrapolate_second_order_and_limit_by_edge()
+                Q.extrapolate_second_order_and_limit_by_vertex()
+            else:
+                raise Exception('Unknown order: %s' % str(self._order_))
+
+
+        w_V     = W.vertex_values
+        u_V     = U.vertex_values
+        v_V     = V.vertex_values
+        z_V     = Z.vertex_values
+
+        h_V     = H.vertex_values
+        uh_V    = UH.vertex_values
+        vh_V    = VH.vertex_values
+
+
+        # Update other quantities
+        #protect(self)
+
+        z_V[:]  = w_V - h_V
+        uh_V[:] = u_V * h_V
+        vh_V[:] = v_V * h_V
+
+
+        num_min = num.min(h_V)
+        if num_min < -1.0e-14:
+            #print 'num.min(h_V)', num_min
+            pass
+
+
+        # Compute edge values by interpolation
+        for name in ['xmomentum', 'ymomentum', 'elevation']:
+            Q = self.quantities[name]
+            Q.interpolate_from_vertices_to_edges()
+
+
+
+
+    def distribute_to_vertices_and_edges_h(self):
+        """Distribution from centroids to edges specific to the SWW eqn.
+
+        It will ensure that h (w-z) is always non-negative even in the
+        presence of steep bed-slopes by taking a weighted average between shallow
+        and deep cases.
+
+        In addition, all conserved quantities get distributed as per either a
+        constant (order==1) or a piecewise linear function (order==2).
+        
+
+        Precondition:
+        All conserved quantities defined at centroids and bed elevation defined at
+        edges.
+        
+        Postcondition
+        Evolved quantities defined at vertices and edges
+        """
+
+
+        #Shortcuts
+        W  = self.quantities['stage']
+        UH = self.quantities['xmomentum']
+        VH = self.quantities['ymomentum']
+        H  = self.quantities['height']
+        Z  = self.quantities['elevation']
+        U  = self.quantities['xvelocity']
+        V  = self.quantities['yvelocity']
+
         #Arrays   
         w_C   = W.centroid_values    

trunk/anuga_work/shallow_water_balanced_steve/swb_domain_ext.c

@@ -612 +612 @@
     k3 = 3*k;  // base index
     
-    // Get stage
+    // Get stage centroid values
     w0 = ((double*) w -> data)[k3 + 0];
     w1 = ((double*) w -> data)[k3 + 1];