Changeset 4815

Timestamp:

Nov 14, 2007, 5:43:02 PM (17 years ago)

Author:

ole

Message:

Work towards Froude_number adjusted, balanced limiters. This is currently
only in place for tight_slope_limiters == 1.

Location:

anuga_core/source/anuga

Files:

• abstract_2d_finite_volumes/quantity.py (modified) (4 diffs)
• config.py (modified) (1 diff)
• shallow_water/shallow_water_ext.c (modified) (6 diffs)
• shallow_water/test_shallow_water_domain.py (modified) (16 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/quantity.py

@@ -24 +24 @@
 from anuga.fit_interpolate.fit import fit_to_mesh
 from anuga.config import points_file_block_line_size as default_block_line_size
+from anuga.config import epsilon
 
 class Quantity:
@@ -109 +110 @@
         (except for a filename or points, attributes (for now))
         - see set_values for details
-
-        Note that if two quantitites q1 and q2 are multiplied,
-        vertex values are multiplied entry by entry
-        while centroid and edge values are re-interpolated.
-        Hence they won't be the product of centroid or edge values
-        from q1 and q2.
-        """
-
-        Q = Quantity(self.domain)
-        Q.set_values(other)
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:    
+            Q = Quantity(self.domain)
+            Q.set_values(other)
 
         result = Quantity(self.domain)
-        result.set_values(self.vertex_values * Q.vertex_values)
+
+        # The product of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values (which calls interpolate). Otherwise
+        # edge and centroid values wouldn't be products from q1 and q2
+        result.vertex_values = self.vertex_values * Q.vertex_values
+        result.edge_values = self.edge_values * Q.edge_values
+        result.centroid_values = self.centroid_values * Q.centroid_values
+        
         return result
 
@@ -129 +135 @@
         return self * other
 
+    def __div__(self, other):
+        """Divide self with anything that could populate a quantity
+
+        E.g other can be a constant, an array, a function, another quantity
+        (except for a filename or points, attributes (for now))
+        - see set_values for details
+
+        Zero division is dealt with by adding an epsilon to the divisore
+        FIXME (Ole): Replace this with native INF once we migrate to NumPy
+        """
+
+        if isinstance(other, Quantity):
+            Q = other
+        else:    
+            Q = Quantity(self.domain)
+            Q.set_values(other)
+
+        result = Quantity(self.domain)
+
+        # The quotient of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values (which calls interpolate). Otherwise
+        # edge and centroid values wouldn't be quotient of q1 and q2
+        result.vertex_values = self.vertex_values/(Q.vertex_values + epsilon)
+        result.edge_values = self.edge_values/(Q.edge_values + epsilon)
+        result.centroid_values = self.centroid_values/(Q.centroid_values + epsilon)
+
+        return result
+
+    def __rdiv__(self, other):
+        """Handle cases like 3/Q, where Q is an instance of class Quantity
+        """
+        return self / other
+
     def __pow__(self, other):
         """Raise quantity to (numerical) power
@@ -140 +180 @@
         """
 
+        if isinstance(other, Quantity):
+            Q = other
+        else:    
+            Q = Quantity(self.domain)
+            Q.set_values(other)
+
         result = Quantity(self.domain)
-        result.set_values(self.vertex_values**other)
+
+        # The power of vertex_values, edge_values and centroid_values
+        # are calculated and assigned directly without using
+        # set_values (which calls interpolate). Otherwise
+        # edge and centroid values wouldn't be correct
+        result.vertex_values = self.vertex_values ** other
+        result.edge_values = self.edge_values ** other
+        result.centroid_values = self.centroid_values ** other
+
         return result
 
-
+    #def __sqrt__(self, other):
+    #    """Define in terms of x**0.5
+    #    """
+    #    pass
+        
 
     def interpolate(self):

anuga_core/source/anuga/config.py

@@ -153 +153 @@
 
 
+maximum_froude_number = 100.0 # To be used in limiters.
+
 minimum_storable_height = 1.0e-5 # Water depth below which it is *stored* as 0
 

anuga_core/source/anuga/shallow_water/shallow_water_ext.c

@@ -141 +141 @@
 
 
-void adjust_froude_number(double *uh,
-                          double h, 
-                          double g) {
+double compute_froude_number(double uh,
+                             double h, 
+                             double g,
+                             double epsilon) {
                           
-  // Adjust momentum if Froude number is excessive
-  double max_froude_number = 20.0;                        
+  // Compute Froude number; v/sqrt(gh)
+  
   double froude_number;
   
   //Compute Froude number (stability diagnostics)
-  froude_number = *uh/sqrt(g*h)/h;
-
-  if (froude_number > max_froude_number) {
-    printf("---------------------------------------------\n");
-    printf("froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h);
-    
-    *uh = *uh/fabs(*uh) * max_froude_number * sqrt(g*h)*h;
-    
-    froude_number = *uh/sqrt(g*h)/h;    
-    printf("Adjusted froude_number=%f (uh=%f, h=%f)\n", froude_number, *uh, h);
-    printf("---------------------------------------------\n");    
-  }
+  if (h > epsilon) { 
+    froude_number = uh/sqrt(g*h)/h;
+  } else {
+    froude_number = 0.0;
+    // FIXME (Ole): What should it be when dry??
+  }
+  
+  return froude_number;
 }
 
@@ -176 +173 @@
   double u;
 
-  //adjust_froude_number(uh, *h, 9.81); // Highly experimental and 
-                                        // probably unneccessary
   
   if (*h < epsilon) {
@@ -484 +479 @@
                               double alpha_balance) {
 
-  int k, k3, i;
+  int k, k3, i, excessive_froude_number=0;
   double dz, hmin, alpha, h_diff, hc_k;
-  double epsilon = 1.0e-6; // Temporary measure
-  double hv[3]; // Depths at vertices
+  double epsilon = 1.0e-6; // FIXME: Temporary measure
+  double g = 9.81; // FIXME: Temporary measure
+  double hv[3], h; // Depths at vertices
+  double Fx, Fy; // Froude numbers
 
   // Compute linear combination between w-limited stages and
@@ -511 +508 @@
     }
 
-    // Calculate depth at vertices
+    // Calculate depth at vertices (possibly negative here!)
     hv[0] = wv[k3] -   zv[k3];
     hv[1] = wv[k3+1] - zv[k3+1];
@@ -606 +603 @@
     //   Momentum is balanced between constant and limited
 
-    if (alpha < 1) {
+
+    if (alpha < 1) {      
       for (i=0; i<3; i++) {
-      
         // FIXME (Ole): Simplify when (if) hvbar gets retired       
         if (beta_h > epsilon) {   
@@ -622 +619 @@
         xmomv[k3+i] = (1-alpha)*xmomc[k] + alpha*xmomv[k3+i];
         ymomv[k3+i] = (1-alpha)*ymomc[k] + alpha*ymomv[k3+i];
+      }
+    }
+        
+
+        
+    if (tight_slope_limiters == 1) {                    
+    
+      // Ensure that the Froude number is kept realistic at vertices
+      // FIXME (Ole): I think it could be used to adjust alpha down
+      // whenever Fr is too large. Possible make sure it doesn't deviate 
+      // too much from the value at the centroid. I like this idea!
+      
+      // FIXME (Ole): currently only used with tights_SL
+
+      excessive_froude_number=0;    
+      for (i=0; i<3; i++) {    
+      
+        // Recalculate depth at vertex i
+        h = wv[k3+i] - zv[k3+i];
+        
+        Fx = compute_froude_number(xmomv[k3+i], h, g, epsilon);
+        Fy = compute_froude_number(ymomv[k3+i], h, g, epsilon);
+        
+        if ( (fabs(Fx) > 100.0) || (fabs(Fy) > 100.0)) {
+          // FIXME: use max_froude - or base it on centroid value of F
+          // printf("Excessive Froude number detected: %f or %f\n", Fx, Fy);
+          excessive_froude_number=1;          
+        }
+      }
+     
+      if (excessive_froude_number) { 
+    
+        // printf("Adjusting momentum to first order.\n"); 
+        // Go back to first order (alpha = 0) for this triangle
+        for (i=0; i<3; i++) {          
+          xmomv[k3+i] = xmomc[k];
+          ymomv[k3+i] = ymomc[k];
+          
+          if (beta_h > epsilon) {         
+            wv[k3+i] = zv[k3+i] + hvbar[k3+i];
+          } else {
+            wv[k3+i] = zv[k3+i] + hc_k; 
+          }
+        }       
       }
     }

anuga_core/source/anuga/shallow_water/test_shallow_water_domain.py

@@ -33 +33 @@
     def __call__(self, x, y):
         from Numeric import zeros, Float
-        from math import sqrt
 
         N = len(x)
@@ -176 +175 @@
     """
 
-    from math import sqrt, exp, cos, pi
+    from math import exp, cos, pi
 
     x = array(x)
@@ -206 +205 @@
     """Rotating field
     """
-    from math import sqrt, atan, pi
+    from math import atan, pi
 
     x = array(x)
@@ -366 +365 @@
         assert allclose(max_speed, 1.31414103233)
 
-    def test_flux_computation(self):    
+    def test_flux_computation(self):    
         """test_flux_computation - test flux calculation (actual C implementation)
-        This one tests the constant case where only the pressure term contributes to each edge and cancels out 
-        once the total flux has been summed up.
-        """
-                
+        This one tests the constant case where only the pressure term contributes to each edge and cancels out 
+        once the total flux has been summed up.
+        """
+                
         a = [0.0, 0.0]
         b = [0.0, 2.0]
@@ -386 +385 @@
         domain.check_integrity()
 
-        # The constant case             
-        domain.set_quantity('elevation', -1)
-        domain.set_quantity('stage', 1) 
-        
-        domain.compute_fluxes()
-        assert allclose(domain.get_quantity('stage').explicit_update[1], 0) # Central triangle
-        
-
-        # The more general case                 
+        # The constant case             
+        domain.set_quantity('elevation', -1)
+        domain.set_quantity('stage', 1) 
+        
+        domain.compute_fluxes()
+        assert allclose(domain.get_quantity('stage').explicit_update[1], 0) # Central triangle
+        
+
+        # The more general case                 
         def surface(x,y):
             return -x/2                    
-        
-        domain.set_quantity('elevation', -10)
-        domain.set_quantity('stage', surface)   
-        domain.set_quantity('xmomentum', 1)             
-        
-        domain.compute_fluxes()
-        
-        #print domain.get_quantity('stage').explicit_update
-        # FIXME (Ole): TODO the general case
-        #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
-                
-        
-                
+        
+        domain.set_quantity('elevation', -10)
+        domain.set_quantity('stage', surface)   
+        domain.set_quantity('xmomentum', 1)             
+        
+        domain.compute_fluxes()
+        
+        #print domain.get_quantity('stage').explicit_update
+        # FIXME (Ole): TODO the general case
+        #assert allclose(domain.get_quantity('stage').explicit_update[1], ........??)
+                
+        
+                
     def test_sw_domain_simple(self):
         a = [0.0, 0.0]
@@ -1751 +1750 @@
     def test_constant_wind_stress(self):
         from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin, sqrt
+        from math import pi, cos, sin
 
         a = [0.0, 0.0]
@@ -1803 +1802 @@
     def test_variable_wind_stress(self):
         from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin, sqrt
+        from math import pi, cos, sin
 
         a = [0.0, 0.0]
@@ -1872 +1871 @@
     def test_windfield_from_file(self):
         from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin, sqrt
+        from math import pi, cos, sin
         from anuga.config import time_format
         from anuga.abstract_2d_finite_volumes.util import file_function
@@ -1982 +1981 @@
     def test_windfield_from_file_seconds(self):
         from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin, sqrt
+        from math import pi, cos, sin
         from anuga.config import time_format
         from anuga.abstract_2d_finite_volumes.util import file_function
@@ -2096 +2095 @@
 
         from anuga.config import rho_a, rho_w, eta_w
-        from math import pi, cos, sin, sqrt
+        from math import pi, cos, sin
 
         a = [0.0, 0.0]
@@ -2586 +2585 @@
     def test_balance_deep_and_shallow(self):
         """Test that balanced limiters preserve conserved quantites.
+        This test is using old depth based balanced limiters
         """
         import copy
@@ -2601 +2601 @@
         elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
 
-        mesh = Domain(points, elements)
-        mesh.check_integrity()
+        domain = Domain(points, elements)
+        domain.check_integrity()
 
         #Create a deliberate overshoot
-        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
-        mesh.set_quantity('elevation', 0) #Flat bed
-        stage = mesh.quantities['stage']
+        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        domain.set_quantity('elevation', 0) #Flat bed
+        stage = domain.quantities['stage']
 
         ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
 
         #Limit
-        mesh.distribute_to_vertices_and_edges()
+        domain.tight_slope_limiters = 0                
+        domain.distribute_to_vertices_and_edges()
 
         #Assert that quantities are conserved
         from Numeric import sum
-        for k in range(len(mesh)):
+        for k in range(len(domain)):
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
@@ -2623 +2624 @@
         #Now try with a non-flat bed - closely hugging initial stage in places
         #This will create alphas in the range [0, 0.478260, 1]
-        mesh.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
-        mesh.set_quantity('elevation', [[0,0,0],
+        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        domain.set_quantity('elevation', [[0,0,0],
                                         [1.8,1.9,5.9],
                                         [4.6,0,0],
                                         [0,2,4]])
-        stage = mesh.quantities['stage']
+        stage = domain.quantities['stage']
 
         ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
@@ -2634 +2635 @@
 
         #Limit
-        mesh.distribute_to_vertices_and_edges()
+        domain.tight_slope_limiters = 0        
+        domain.distribute_to_vertices_and_edges()
 
 
         #Assert that all vertex quantities have changed
-        for k in range(len(mesh)):
+        for k in range(len(domain)):
             #print ref_vertex_values[k,:], stage.vertex_values[k,:]
             assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
         #and assert that quantities are still conserved
         from Numeric import sum
-        for k in range(len(mesh)):
+        for k in range(len(domain)):
+            assert allclose (ref_centroid_values[k],
+                             sum(stage.vertex_values[k,:])/3)
+
+
+        # Check actual results
+        assert allclose (stage.vertex_values,
+                         [[2,2,2],
+                          [1.93333333, 2.03333333, 6.03333333],
+                          [6.93333333, 4.53333333, 4.53333333],
+                          [5.33333333, 5.33333333, 5.33333333]])
+
+
+    def test_balance_deep_and_shallow_tight_SL(self):
+        """Test that balanced limiters preserve conserved quantites.
+        This test is using Tight Slope Limiters
+        """
+        import copy
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+
+        #bac, bce, ecf, dbe
+        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        domain = Domain(points, elements)
+        domain.check_integrity()
+
+        #Create a deliberate overshoot
+        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        domain.set_quantity('elevation', 0) #Flat bed
+        stage = domain.quantities['stage']
+
+        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
+
+        #Limit
+        domain.tight_slope_limiters = 1                
+        domain.distribute_to_vertices_and_edges()
+
+        #Assert that quantities are conserved
+        from Numeric import sum
+        for k in range(len(domain)):
+            assert allclose (ref_centroid_values[k],
+                             sum(stage.vertex_values[k,:])/3)
+
+
+        #Now try with a non-flat bed - closely hugging initial stage in places
+        #This will create alphas in the range [0, 0.478260, 1]
+        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        domain.set_quantity('elevation', [[0,0,0],
+                                        [1.8,1.9,5.9],
+                                        [4.6,0,0],
+                                        [0,2,4]])
+        stage = domain.quantities['stage']
+
+        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
+        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
+
+        #Limit
+        domain.tight_slope_limiters = 1        
+        domain.distribute_to_vertices_and_edges()
+
+
+        #Assert that all vertex quantities have changed
+        for k in range(len(domain)):
+            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
+            assert not allclose (ref_vertex_values[k,:], stage.vertex_values[k,:])
+        #and assert that quantities are still conserved
+        from Numeric import sum
+        for k in range(len(domain)):
             assert allclose (ref_centroid_values[k],
                              sum(stage.vertex_values[k,:])/3)
@@ -2657 +2734 @@
         #                  [5.33333333, 5.33333333, 5.33333333]])
 
+
+
+    def test_balance_deep_and_shallow_Froude(self):
+        """Test that balanced limiters preserve conserved quantites -
+        and also that excessive Froude numbers are dealt with.
+        This test is using tight slope limiters.
+        """
+        import copy
+        from Numeric import sqrt, absolute
+
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+
+        # bac, bce, ecf, dbe
+        elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ]
+
+        domain = Domain(points, elements)
+        domain.check_integrity()
+
+        # Create non-flat bed - closely hugging initial stage in places
+        # This will create alphas in the range [0, 0.478260, 1]
+        domain.set_quantity('stage', [[3,0,3], [2,2,6], [5,3,8], [8,3,5]])
+        domain.set_quantity('elevation', [[0,0,0],
+                                        [1.8,1.999,5.999],
+                                        [4.6,0,0],
+                                        [0,2,4]])
+
+        # Create small momenta, that nonetheless will generate large speeds
+        # due to shallow depth at isolated vertices
+        domain.set_quantity('xmomentum', -0.0058)
+        domain.set_quantity('ymomentum', 0.0890)
+
+
+
+        
+        stage = domain.quantities['stage']
+        elevation = domain.quantities['elevation']
+        xmomentum = domain.quantities['xmomentum']
+        ymomentum = domain.quantities['ymomentum']        
+
+        # Setup triangle #1 to mimick real Froude explosion observed
+        # in the Onslow example 13 Nov 2007.
+
+        stage.vertex_values[1,:] = [1.6385, 1.6361, 1.2953]
+        elevation.vertex_values[1,:] = [1.6375, 1.6336, 0.4647]        
+        xmomentum.vertex_values[1,:] = [-0.0058, -0.0050, -0.0066]
+        ymomentum.vertex_values[1,:] = [0.0890, 0.0890, 0.0890]
+
+        xmomentum.interpolate()
+        ymomentum.interpolate()        
+        stage.interpolate()       
+        elevation.interpolate()
+
+        # Verify interpolation
+        assert allclose(stage.centroid_values[1], 1.5233)
+        assert allclose(elevation.centroid_values[1], 1.2452667)
+        assert allclose(xmomentum.centroid_values[1], -0.0058)
+        assert allclose(ymomentum.centroid_values[1], 0.089)
+
+        # Derived quantities
+        depth = stage-elevation
+        u = xmomentum/depth
+        v = ymomentum/depth
+
+        denom = (depth*g)**0.5 
+        Fx = u/denom
+        Fy = v/denom
+        
+   
+        # Verify against Onslow example (14 Nov 2007)
+        assert allclose(depth.centroid_values[1], 0.278033)
+        assert allclose(u.centroid_values[1], -0.0208608)
+        assert allclose(v.centroid_values[1], 0.3201055)
+
+        assert allclose(denom.centroid_values[1],
+                        sqrt(depth.centroid_values[1]*g))
+
+        assert allclose(u.centroid_values[1]/denom.centroid_values[1],
+                        -0.012637746977)
+        assert allclose(Fx.centroid_values[1],
+                        u.centroid_values[1]/denom.centroid_values[1])
+
+        # Check that Froude numbers are small at centroids.
+        assert allclose(Fx.centroid_values[1], -0.012637746977)
+        assert allclose(Fy.centroid_values[1], 0.193924048435)
+
+
+        # But Froude numbers are huge at some vertices and edges
+        assert allclose(Fx.vertex_values[1,:], [-5.85888475e+01,
+                                                -1.27775313e+01,
+                                                -2.78511420e-03])
+
+        assert allclose(Fx.edge_values[1,:], [-6.89150773e-03,
+                                              -7.38672488e-03,
+                                              -2.35626238e+01])
+
+        assert allclose(Fy.vertex_values[1,:], [8.99035764e+02,
+                                                2.27440057e+02,
+                                                3.75568430e-02])
+
+        assert allclose(Fy.edge_values[1,:], [1.05748998e-01,
+                                              1.06035244e-01,
+                                              3.88346947e+02])
+        
+        
+        # The task is now to arrange the limiters such that Froude numbers
+        # remain under control whil at the same time obeying the conservation
+        # laws.
+
+        
+        ref_centroid_values = copy.copy(stage.centroid_values[:]) #Copy
+        ref_vertex_values = copy.copy(stage.vertex_values[:]) #Copy
+
+        # Limit (and invoke balance_deep_and_shallow)
+        domain.tight_slope_limiters = 1
+        domain.distribute_to_vertices_and_edges()
+
+        # Redo derived quantities
+        depth = stage-elevation
+        u = xmomentum/depth
+        v = ymomentum/depth
+
+        denom = (depth*g)**0.5 
+        Fx = u/denom
+        Fy = v/denom
+
+
+        # Assert that Froude numbers are less than max value (TBA)
+        # at vertices, edges and centroids.
+        from anuga.config import maximum_froude_number
+        assert alltrue(absolute(Fx.vertex_values[1,:]) < maximum_froude_number)
+        assert alltrue(absolute(Fy.vertex_values[1,:]) < maximum_froude_number)
+
+
+        # Assert that all vertex quantities have changed
+        for k in range(len(domain)):
+            #print ref_vertex_values[k,:], stage.vertex_values[k,:]
+            assert not allclose (ref_vertex_values[k,:],
+                                 stage.vertex_values[k,:])
+            
+        # Assert that quantities are still conserved
+        from Numeric import sum
+        for k in range(len(domain)):
+            assert allclose (ref_centroid_values[k],
+                             sum(stage.vertex_values[k,:])/3)
+
+
+        
+        return
+    
+        qwidth = 12
+        for k in [1]: #range(len(domain)):
+            print 'Triangle %d (C, V, E)' %k
+            
+            print 'stage'.ljust(qwidth), stage.centroid_values[k],\
+                  stage.vertex_values[k,:], stage.edge_values[k,:]
+            print 'elevation'.ljust(qwidth), elevation.centroid_values[k],\
+                  elevation.vertex_values[k,:], elevation.edge_values[k,:]
+            print 'depth'.ljust(qwidth), depth.centroid_values[k],\
+                  depth.vertex_values[k,:], depth.edge_values[k,:]
+            print 'xmomentum'.ljust(qwidth), xmomentum.centroid_values[k],\
+                  xmomentum.vertex_values[k,:], xmomentum.edge_values[k,:]
+            print 'ymomentum'.ljust(qwidth), ymomentum.centroid_values[k],\
+                  ymomentum.vertex_values[k,:], ymomentum.edge_values[k,:]
+            print 'u'.ljust(qwidth),u.centroid_values[k],\
+                  u.vertex_values[k,:], u.edge_values[k,:]
+            print 'v'.ljust(qwidth), v.centroid_values[k],\
+                  v.vertex_values[k,:], v.edge_values[k,:]
+            print 'Fx'.ljust(qwidth), Fx.centroid_values[k],\
+                  Fx.vertex_values[k,:], Fx.edge_values[k,:]
+            print 'Fy'.ljust(qwidth), Fy.centroid_values[k],\
+                  Fy.vertex_values[k,:], Fy.edge_values[k,:]
+            
+            
+        
 
 
@@ -5306 +5564 @@
 
     suite = unittest.makeSuite(Test_Shallow_Water,'test')
+
+    #suite = unittest.makeSuite(Test_Shallow_Water,'test_balance_deep_and_shallow_Froude')
+    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_fitting_using_shallow_water_domain')    
     #suite = unittest.makeSuite(Test_Shallow_Water,'test_tight_slope_limiters')