Changeset 7573 for anuga_core/source/anuga/shallow_water_balanced/swb_domain.py

Timestamp:

Nov 27, 2009, 9:31:34 AM (14 years ago)

Author:

steve

Message:

Committing a version of shallow_water_balanced which passes it unit test
using a version of edge limiting which doesn't limit boundary edges. THis
is useful to allow linear functions to be reconstructed.

Had to play with the transmissive BC and use centroid values which is
set via domain set_centroid_transmissive_bc

File:

• anuga_core/source/anuga/shallow_water_balanced/swb_domain.py (modified) (6 diffs)

anuga_core/source/anuga/shallow_water_balanced/swb_domain.py

@@ -3 +3 @@
 from anuga.shallow_water.shallow_water_domain import Domain as Sww_domain
 
+from swb_boundary_conditions import Transmissive_boundary
 
 ##############################################################################
@@ -15 +16 @@
 
     ##
-    # @brief Instantiate a shallow water domain.
+    # @brief Instantiate a shallow water balanced domain.
     # @param coordinates
     # @param vertices
@@ -48 +49 @@
                  number_of_full_triangles=None):
 
+        conserved_quantities = [ 'stage', 'xmomentum', 'ymomentum']
+
         evolved_quantities = [ 'stage', 'xmomentum', 'ymomentum', \
                                'height', 'elevation', 'xvelocity', 'yvelocity']
@@ -64 +67 @@
                             use_cache = use_cache,
                             verbose = verbose,
+                            conserved_quantities = conserved_quantities,
                             evolved_quantities = evolved_quantities,
                             other_quantities = other_quantities,
@@ -76 +80 @@
         # set some defaults
         #---------------------
-        self.set_timestepping_method('euler')
-        self.set_default_order(1)
+        self.set_timestepping_method(1)
+        self.set_default_order(2)
         self.set_new_mannings_function(True)
-        self.set_use_edge_limiter(True)
-
-
-
+        self.set_centroid_transmissive_bc(True)
+
+    ##
+    # @brief Run integrity checks on shallow water balanced domain.
+    def check_integrity(self):
+        Sww_domain.check_integrity(self)
+
+        #Check that the evolved quantities are correct (order)
+        msg = 'First evolved quantity must be "stage"'
+        assert self.evolved_quantities[0] == 'stage', msg
+        msg = 'Second evolved quantity must be "xmomentum"'
+        assert self.evolved_quantities[1] == 'xmomentum', msg
+        msg = 'Third evolved quantity must be "ymomentum"'
+        assert self.evolved_quantities[2] == 'ymomentum', msg
+        msg = 'Fourth evolved quantity must be "height"'
+        assert self.evolved_quantities[3] == 'height', msg
+        msg = 'Fifth evolved quantity must be "elevation"'
+        assert self.evolved_quantities[4] == 'elevation', msg
+        msg = 'Sixth evolved quantity must be "xvelocity"'
+        assert self.evolved_quantities[5] == 'xvelocity', msg        
+        msg = 'Seventh evolved quantity must be "yvelocity"'
+        assert self.evolved_quantities[6] == 'yvelocity', msg        
+
+        msg = 'First other quantity must be "friction"'
+        assert self.other_quantities[0] == 'friction', msg
+
+    ##
+    # @brief 
+    def compute_fluxes(self):
+        #Call correct module function (either from this module or C-extension)
+        compute_fluxes(self)
+
+    ##
+    # @brief 
+    def distribute_to_vertices_and_edges(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
+        Stage  = self.quantities['stage']
+        Xmom   = self.quantities['xmomentum']
+        Ymom   = self.quantities['ymomentum']
+        Elev   = self.quantities['elevation']
+        Height = self.quantities['height']
+        Xvel   = self.quantities['xvelocity']
+        Yvel   = self.quantities['yvelocity']
+
+        #Arrays   
+        w_C   = Stage.centroid_values    
+        uh_C  = Xmom.centroid_values
+        vh_C  = Ymom.centroid_values    
+        z_C   = Elev.centroid_values
+        h_C   = Height.centroid_values
+        u_C   = Xvel.centroid_values
+        v_C   = Yvel.centroid_values
+
+        w_C[:] = num.maximum(w_C, z_C)
+        
+        h_C[:] = w_C - z_C
+
+
+        assert num.min(h_C) >= 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-15)        
+        
+        u_C[:]  = uh_C/h_C
+        v_C[:]  = vh_C/h_C
+        
+        for name in [ 'height', '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()
+            else:
+                raise 'Unknown order'
+
+
+        w_E     = Stage.edge_values
+        uh_E    = Xmom.edge_values
+        vh_E    = Ymom.edge_values      
+        z_E     = Elev.edge_values      
+        h_E     = Height.edge_values
+        u_E     = Xvel.edge_values
+        v_E     = Yvel.edge_values              
+
+
+        w_E[:]   = z_E + h_E
+
+        #num.putmask(u_E, h_temp <= 0.0, 0.0)
+        #num.putmask(v_E, h_temp <= 0.0, 0.0)
+        #num.putmask(w_E, h_temp <= 0.0, z_E+h_E)
+        #num.putmask(h_E, h_E <= 0.0, 0.0)
+        
+        uh_E[:] = u_E * h_E
+        vh_E[:] = v_E * h_E
+
+        """
+        print '=========================================================='
+        print 'Time ', self.get_time()
+        print h_E
+        print uh_E
+        print vh_E
+        """
+        
+        # Compute vertex values by interpolation
+        for name in self.evolved_quantities:
+            Q = self.quantities[name]
+            Q.interpolate_from_edges_to_vertices()
+
+
+        w_V     = Stage.vertex_values
+        uh_V    = Xmom.vertex_values
+        vh_V    = Ymom.vertex_values    
+        z_V     = Elev.vertex_values    
+        h_V     = Height.vertex_values
+        u_V     = Xvel.vertex_values
+        v_V     = Yvel.vertex_values            
+
+
+        #w_V[:]    = z_V + h_V
+
+        #num.putmask(u_V, h_V <= 0.0, 0.0)
+        #num.putmask(v_V, h_V <= 0.0, 0.0)
+        #num.putmask(w_V, h_V <= 0.0, z_V)        
+        #num.putmask(h_V, h_V <= 0.0, 0.0)
+        
+        uh_V[:] = u_V * h_V
+        vh_V[:] = v_V * h_V
+
+
+    ##
+    # @brief Code to let us use old shallow water domain BCs
     def conserved_values_to_evolved_values(self, q_cons, q_evol):
         """Mapping between conserved quantities and the evolved quantities.
@@ -111 +262 @@
         hc = wc - be
 
-        if hc < 0.0:
+        if hc <= 0.0:
             hc = 0.0
             uc = 0.0