Changeset 6689 for branches/numpy/anuga/shallow_water/shallow_water_domain.py

Timestamp:

Apr 1, 2009, 3:19:07 PM (15 years ago)

Author:

rwilson

Message:

Back-merge from Numeric trunk to numpy branch.

File:

• branches/numpy/anuga/shallow_water/shallow_water_domain.py (modified) (7 diffs)

branches/numpy/anuga/shallow_water/shallow_water_domain.py

@@ -812 +812 @@
 
         return msg
-
+       
+        
+
+    def compute_boundary_flows(self):
+        """Compute boundary flows at current timestep.
+                        
+        Quantities computed are:
+           Total inflow across boundary
+           Total outflow across boundary
+           Flow across each tagged boundary segment
+        """
+                
+        # Run through boundary array and compute for each segment
+        # the normal momentum ((uh, vh) dot normal) times segment length.
+        # Based on sign accumulate this into boundary_inflow and boundary_outflow.
+                        
+        # Compute flows along boundary
+        
+        uh = self.get_quantity('xmomentum').get_values(location='edges')
+        vh = self.get_quantity('ymomentum').get_values(location='edges')        
+        
+        # Loop through edges that lie on the boundary and calculate 
+        # flows
+        boundary_flows = {}
+        total_boundary_inflow = 0.0
+        total_boundary_outflow = 0.0
+        for vol_id, edge_id in self.boundary:
+            # Compute normal flow across edge. Since normal vector points
+            # away from triangle, a positive sign means that water 
+            # flows *out* from this triangle.
+             
+            momentum = [uh[vol_id, edge_id], vh[vol_id, edge_id]]
+            normal = self.mesh.get_normal(vol_id, edge_id)
+            length = self.mesh.get_edgelength(vol_id, edge_id)            
+            normal_flow = num.dot(momentum, normal)*length
+            
+            # Reverse sign so that + is taken to mean inflow
+            # and - means outflow. This is more intuitive.
+            edge_flow = -normal_flow
+            
+            # Tally up inflows and outflows separately
+            if edge_flow > 0:
+                # Flow is inflow      
+                total_boundary_inflow += edge_flow                                  
+            else:
+                # Flow is outflow
+                total_boundary_outflow += edge_flow    
+
+            # Tally up flows by boundary tag
+            tag = self.boundary[(vol_id, edge_id)]
+            
+            if tag not in boundary_flows:
+                boundary_flows[tag] = 0.0
+            boundary_flows[tag] += edge_flow
+            
+                
+        return boundary_flows, total_boundary_inflow, total_boundary_outflow
+        
+
+    def compute_forcing_flows(self):
+        """Compute flows in and out of domain due to forcing terms.
+                        
+        Quantities computed are:
+                
+        
+           Total inflow through forcing terms
+           Total outflow through forcing terms
+           Current total volume in domain        
+
+        """
+
+        #FIXME(Ole): We need to separate what part of explicit_update was 
+        # due to the normal flux calculations and what is due to forcing terms.
+        
+        pass
+                        
+        
+    def compute_total_volume(self):
+        """Compute total volume (m^3) of water in entire domain
+        """
+        
+        area = self.mesh.get_areas()
+        volume = 0.0
+        
+        stage = self.get_quantity('stage').get_values(location='centroids')
+        elevation = self.get_quantity('elevation').get_values(location='centroids')        
+        depth = stage-elevation
+        
+        #print 'z', elevation
+        #print 'w', stage
+        #print 'h', depth
+        return num.sum(depth*area)
+        
+        
+    def volumetric_balance_statistics(self):                
+        """Create volumetric balance report suitable for printing or logging.
+        """
+        
+        boundary_flows, total_boundary_inflow, total_boundary_outflow = self.compute_boundary_flows() 
+        
+        s = '---------------------------\n'        
+        s += 'Volumetric balance report:\n'
+        s += '--------------------------\n'
+        s += 'Total boundary inflow [m^3/s]: %.2f\n' % total_boundary_inflow
+        s += 'Total boundary outflow [m^3/s]: %.2f\n' % total_boundary_outflow        
+        s += 'Net boundary flow by tags [m^3/s]\n'
+        for tag in boundary_flows:
+            s += '    %s [m^3/s]: %.2f\n' % (tag, boundary_flows[tag])
+        
+        s += 'Total net boundary flow [m^3/s]: %.2f\n' % (total_boundary_inflow + total_boundary_outflow) 
+        s += 'Total volume in domain [m^3]: %.2f\n' % self.compute_total_volume()
+        
+        # The go through explicit forcing update and record the rate of change for stage and 
+        # record into forcing_inflow and forcing_outflow. Finally compute integral 
+        # of depth to obtain total volume of domain.
+        
+        # FIXME(Ole): This part is not yet done.                
+        
+        return s        
+           
 ################################################################################
 # End of class Shallow Water Domain
@@ -1386 +1505 @@
         domain: pointer to shallow water domain for which the boundary applies
         mean_stage: The mean water level which will be added to stage derived
-                    from the sww file
+                    from the boundary condition
         time_thinning: Will set how many time steps from the sww file read in
                        will be interpolated to the boundary. For example if
@@ -1402 +1521 @@
                           that available in the underlying data.
 
+                          Note that mean_stage will also be added to this.
+                                                
         use_cache:
         verbose:
@@ -1574 +1695 @@
                 xmom_update[k] += S*uh[k]
                 ymom_update[k] += S*vh[k]
+
+def depth_dependent_friction(domain, default_friction,
+                             surface_roughness_data,
+                             verbose=False):
+    """Returns an array of friction values for each wet element adjusted for depth.
+
+    Inputs:
+        domain - computational domain object
+        default_friction - depth independent bottom friction
+        surface_roughness_data - N x 5 array of n0, d1, n1, d2, n2 values for each
+        friction region.
+
+    Outputs:
+        wet_friction - Array that can be used directly to update friction as follows:
+                       domain.set_quantity('friction', wet_friction)
+
+        
+        
+    """
+
+    import Numeric as num
+    
+    # Create a temp array to store updated depth dependent friction for wet elements
+    # EHR this is outwardly inneficient but not obvious how to avoid recreating each call??????
+    N=len(domain)
+    wet_friction    = num.zeros(N, num.Float)
+    wet_friction[:] = default_n0   # Initially assign default_n0 to all array so sure have no zeros values
+    
+    
+    depth = domain.create_quantity_from_expression('stage - elevation')  # create depth instance for this timestep
+    # Recompute depth as vector  
+    d = depth.get_values(location='centroids')
+ 
+    # rebuild the 'friction' values adjusted for depth at this instant
+    for i in domain.get_wet_elements():                                  # loop for each wet element in domain
+        
+        # Get roughness data for each element
+        n0 = float(surface_roughness_data[i,0])
+        d1 = float(surface_roughness_data[i,1])
+        n1 = float(surface_roughness_data[i,2])
+        d2 = float(surface_roughness_data[i,3])
+        n2 = float(surface_roughness_data[i,4])
+        
+        
+        # Recompute friction values from depth for this element 
+               
+        if d[i]   <= d1: 
+            depth_dependent_friction = n1
+        elif d[i] >= d2:
+            depth_dependent_friction = n2
+        else:
+            depth_dependent_friction = n1+((n2-n1)/(d2-d1))*(d[i]-d1)
+            
+        # check sanity of result
+        if (depth_dependent_friction  < 0.010 or depth_dependent_friction > 9999.0) :
+            print model_data.basename+' >>>> WARNING: computed depth_dependent friction out of range ddf,n1,n2 ', depth_dependent_friction, n1,n2
+        
+        # update depth dependent friction  for that wet element
+        wet_friction[i] = depth_dependent_friction
+        
+    # EHR add code to show range of 'friction across domain at this instant as sanity check?????????
+    
+    if verbose :
+        nvals=domain.get_quantity('friction').get_values(location='centroids')        # return array of domain nvals
+        n_min=min(nvals)
+        n_max=max(nvals)
+        
+        print "         ++++ calculate_depth_dependent_friction  - Updated friction - range  %7.3f to %7.3f" %(n_min,n_max)
+    
+    return wet_friction
+
 
 ################################################################################
@@ -1901 +2093 @@
                 assert is_inside_polygon(point, bounding_polygon), msg
 
-            # Compute area and check that it is greater than 0
-            self.exchange_area = polygon_area(self.polygon)
-
-            msg = 'Polygon %s in forcing term' % str(self.polygon)
-            msg += ' has area = %f' % self.exchange_area
-            assert self.exchange_area > 0.0
-
         # Pointer to update vector
         self.update = domain.quantities[self.quantity_name].explicit_update
@@ -1931 +2116 @@
         if self.polygon is not None:
             # Inlet is polygon
-            inlet_region = ('polygon=\n%s, area=%f m^2' %
-                            (self.polygon, self.exchange_area))
-
+            inlet_region = 'polygon=%s' % (self.polygon) 
             self.exchange_indices = inside_polygon(points, self.polygon)
 
-        if self.exchange_indices is not None:
+        if self.exchange_indices is None:
+            self.exchange_area = polygon_area(bounding_polygon)
+        else:    
             if len(self.exchange_indices) == 0:
                 msg = 'No triangles have been identified in '
@@ -1942 +2127 @@
                 raise Exception, msg
 
+            # Compute exchange area as the sum of areas of triangles identified
+            # by circle or polygon
+            self.exchange_area = 0.0
+            for i in self.exchange_indices:
+                self.exchange_area += domain.areas[i]
+            
+
+        msg = 'Exchange area in forcing term'
+        msg += ' has area = %f' %self.exchange_area
+        assert self.exchange_area > 0.0            
+            
+                
+
+            
         # Check and store default_rate
         msg = ('Keyword argument default_rate must be either None '