Changeset 7176 for branches/numpy/anuga/shallow_water

Timestamp:

Jun 10, 2009, 5:28:35 PM (16 years ago)

Author:

rwilson

Message:

Back-merge from Numeric to numpy.

Location:

branches/numpy/anuga/shallow_water

Files:

• data_manager.py (modified) (4 diffs)
• shallow_water_domain.py (modified) (1 diff)

branches/numpy/anuga/shallow_water/data_manager.py

@@ -566 +566 @@
                 # Define a zero vector of same size and type as A
                 # for use with momenta
-                null = num.zeros(num.size(A), A.dtype.char) #??#
+                null = num.zeros(num.size(A), A.dtype.char)
 
                 # Get xmomentum where depth exceeds minimum_storable_height
@@ -1433 +1433 @@
     for name in infile.variables:
         var = infile.variables[name]
-        outfile.createVariable(name, var.dtype.char, var.dimensions)    #??#
+        outfile.createVariable(name, var.dtype.char, var.dimensions)
 
     # Copy the static variables
@@ -2193 +2193 @@
     swwfile = basename_in + '.sww'
     demfile = basename_out + '.' + format
-    # Note the use of a .ers extension is optional (write_ermapper_grid will
-    # deal with either option
 
     # Read sww file
@@ -2697 +2695 @@
     # Interpolate using quantity values
     if verbose: print 'Interpolating'
-    interpolated_values = interp.interpolate(q, data_points).flatten()    #????#
+    interpolated_values = interp.interpolate(q, data_points).flatten()
 
     if verbose:

branches/numpy/anuga/shallow_water/shallow_water_domain.py

@@ -1462 +1462 @@
 
 
-##
-# @brief A class for a 'field' boundary.
-# @note Inherits from Boundary.
+class Inflow_boundary(Boundary):
+    """Apply given flow in m^3/s to boundary segment.
+    Depth and momentum is derived using Manning's formula.
+
+    Underlying domain must be specified when boundary is instantiated
+    """
+    
+    # FIXME (Ole): This is work in progress and definitely not finished.
+    # The associated test has been disabled
+
+    def __init__(self, domain=None, rate=0.0):
+        Boundary.__init__(self)
+
+        if domain is None:
+            msg = 'Domain must be specified for '
+            msg += 'Inflow boundary'
+            raise Exception, msg
+
+        self.domain = domain
+        
+        # FIXME(Ole): Allow rate to be time dependent as well
+        self.rate = rate
+        self.tag = None # Placeholder for tag associated with this object.
+
+    def __repr__(self):
+        return 'Inflow_boundary(%s)' %self.domain
+
+    def evaluate(self, vol_id, edge_id):
+        """Apply inflow rate at each edge of this boundary
+        """
+        
+        # First find all segments having the same tag is vol_id, edge_id
+        # This will be done the first time evaluate is called.
+        if self.tag is None:
+            boundary = self.domain.boundary
+            self.tag = boundary[(vol_id, edge_id)]        
+            
+            # Find total length of boundary with this tag
+            length = 0.0
+            for v_id, e_id in boundary:
+                if self.tag == boundary[(v_id, e_id)]:
+                    length += self.domain.mesh.get_edgelength(v_id, e_id)            
+
+            self.length = length
+            self.average_momentum = self.rate/length
+            
+            
+        # Average momentum has now been established across this boundary
+        # Compute momentum in the inward normal direction 
+        
+        inward_normal = -self.domain.mesh.get_normal(vol_id, edge_id)       
+        xmomentum, ymomentum = self.average_momentum * inward_normal
+            
+        # Compute depth based on Manning's formula v = 1/n h^{2/3} sqrt(S)
+        # Where v is velocity, n is manning's coefficient, h is depth and S is the slope into the domain. 
+        # Let mu be the momentum (vh), then this equation becomes: mu = 1/n h^{5/3} sqrt(S) 
+        # from which we can isolate depth to get
+        # h = (mu n/sqrt(S) )^{3/5} 
+        
+        slope = 0 # get gradient for this triangle dot normal
+        
+        # get manning coef from this triangle
+        friction = self.domain.get_quantity('friction').get_values(location='edges', 
+                                                                   indices=[vol_id])[0]
+        mannings_n = friction[edge_id]
+
+        if slope > epsilon and mannings_n > epsilon:
+            depth = pow(self.average_momentum * mannings_n/math.sqrt(slope), 3.0/5) 
+        else:
+            depth = 1.0
+            
+        # Elevation on this edge    
+        
+        z = self.domain.get_quantity('elevation').get_values(location='edges', 
+                                                             indices=[vol_id])[0]
+        elevation = z[edge_id]
+            
+        # Assign conserved quantities and return
+        q = num.array([elevation + depth, xmomentum, ymomentum], num.Float)
+        return q
+
+
+        
+    
+            
+        
 class Field_boundary(Boundary):
     """Set boundary from given field represented in an sww file containing