Changeset 6556 for anuga_core

Timestamp:

Mar 19, 2009, 4:08:50 PM (16 years ago)

Author:

rwilson

Message:

Changes to sww2dem() to minimize memory footprint.

Location:

anuga_core/source/anuga

Files:

• shallow_water/data_manager.py (modified) (6 diffs)
• shallow_water/test_data_manager.py (modified) (1 diff)
• utilities/system_tools.py (modified) (1 diff)
• utilities/test_system_tools.py (modified) (1 diff)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -89 +89 @@
 from anuga.load_mesh.loadASCII import export_mesh_file
 from anuga.utilities.polygon import intersection
-
+from anuga.utilities.system_tools import get_vars_in_expression
+
+
+# Default block size for sww2dem()
+DEFAULT_BLOCK_SIZE = 10000
 
 ######
@@ -2104 +2108 @@
             origin=None,
             datum='WGS84',
-            format='ers'):
+            format='ers',
+            block_size=None):
     """Read SWW file and convert to Digitial Elevation model format
     (.asc or .ers)
@@ -2147 +2152 @@
 
     format can be either 'asc' or 'ers'
+
+    block_size - sets the number of slices along the non-time axis to
+                 process in one block.
     """
 
@@ -2177 +2185 @@
         number_of_decimal_places = 3
 
+    if block_size is None:
+        block_size = DEFAULT_BLOCK_SIZE
+
+    # Read SWW file
     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
@@ -2253 +2263 @@
             if verbose: print '    %s in [%f, %f]' %(name, min(q), max(q))
 
-    # Get quantity and reduce if applicable
-    if verbose: print 'Processing quantity %s' %quantity
-
-    # Turn NetCDF objects into Numeric arrays
-    try:
-        q = fid.variables[quantity][:]
-    except:
-        quantity_dict = {}
-        for name in fid.variables.keys():
-            quantity_dict[name] = fid.variables[name][:]
-        #Convert quantity expression to quantities found in sww file
-        q = apply_expression_to_dictionary(quantity, quantity_dict)
-
-    if len(q.shape) == 2:
-        #q has a time component, must be reduced alongthe temporal dimension
-        if verbose: print 'Reducing quantity %s' %quantity
-        q_reduced = num.zeros(number_of_points, num.Float)
-
-        if timestep is not None:
-            for k in range(number_of_points):
-                q_reduced[k] = q[timestep,k]
-        else:
-            for k in range(number_of_points):
-                q_reduced[k] = reduction( q[:,k] )
-
-        q = q_reduced
-
+    # Get the variables in the supplied expression.
+    # This may throw a SyntaxError exception.
+    var_list = get_vars_in_expression(quantity)
+
+    # Check that we have the required variables in the SWW file.
+    missing_vars = []
+    for name in var_list:
+        try:
+            _ = fid.variables[name]
+        except:
+            missing_vars.append(name)
+    if missing_vars:
+        msg = ("In expression '%s', variables %s are not in the SWW file '%s'"
+               % (quantity, swwfile))
+        raise Exception, msg
+
+    # Create result array and start filling, block by block.
+    result = num.zeros(number_of_points, num.Float)
+
+    for start_slice in xrange(0, number_of_points, block_size):
+        # limit slice size to array end if at last block
+        end_slice = min(start_slice + block_size, number_of_points)
+        
+        # get slices of all required variables
+        q_dict = {}
+        for name in var_list:
+            # check if variable has time axis
+            if len(fid.variables[name].shape) == 2:
+                q_dict[name] = fid.variables[name][:,start_slice:end_slice]
+            else:       # no time axis
+                q_dict[name] = fid.variables[name][start_slice:end_slice]
+
+        # Evaluate expression with quantities found in SWW file
+        res = apply_expression_to_dictionary(quantity, q_dict)
+
+        if len(res.shape) == 2:
+            new_res = num.zeros(res.shape[1], num.Float)
+            for k in xrange(res.shape[1]):
+                new_res[k] = reduction(res[:,k])
+            res = new_res
+
+        result[start_slice:end_slice] = res
+                                    
     #Post condition: Now q has dimension: number_of_points
-    assert len(q.shape) == 1
-    assert q.shape[0] == number_of_points
+    assert len(result.shape) == 1
+    assert result.shape[0] == number_of_points
 
     if verbose:
         print 'Processed values for %s are in [%f, %f]' % \
-              (quantity, min(q), max(q))
+              (quantity, min(result), max(result))
 
     #Create grid and update xll/yll corner and x,y
@@ -2359 +2385 @@
     #Interpolate using quantity values
     if verbose: print 'Interpolating'
-    grid_values = interp.interpolate(q, grid_points).flat
+    grid_values = interp.interpolate(result, grid_points).flat
 
     if verbose:

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -2000 +2000 @@
                 number_of_decimal_places = 9,
                 verbose = self.verbose,
-                format = 'asc')
+                format = 'asc',
+                block_size=2)
 
 

anuga_core/source/anuga/utilities/system_tools.py

@@ -247 +247 @@
         
             
-
-
+##
+# @brief Get list of variable names in an expression string.
+# @param source A string containing a python expression.
+# @return A list of variable name strings.
+# @note Throws SyntaxError exception if not a valid expression.
+def get_vars_in_expression(source):
+    '''Get list of variable names in a python expression.'''
+
+    import compiler
+    from compiler.ast import Node
+
+    ##
+    # @brief Internal recursive function.
+    # @param node An AST parse Node.
+    # @param var_list Input list of variables.
+    # @return An updated list of variables.
+    def get_vars_body(node, var_list=[]):
+        if isinstance(node, Node):
+            if node.__class__.__name__ == 'Name':
+                for child in node.getChildren():
+                    if child not in var_list:
+                        var_list.append(child)
+            if any(isinstance(child, Node) for child in node.getChildren()):
+                for child in node.getChildren():
+                    var_list = get_vars_body(child, var_list)
+
+        return var_list
+
+    return get_vars_body(compiler.parse(source))

anuga_core/source/anuga/utilities/test_system_tools.py

@@ -134 +134 @@
         #print checksum
         
-        
+
+    def test_get_vars_in_expression(self):
+        '''Test the 'get vars from expression' code.'''
+
+        def test_it(source, expected):
+            result = get_vars_in_expression(source)
+            result.sort()
+            expected.sort()
+            msg = ("Source: '%s'\nResult: %s\nExpected: %s"
+                   % (source, str(result), str(expected)))
+            self.failUnlessEqual(result, expected, msg)
+                
+        source = 'fred'
+        expected = ['fred']
+        test_it(source, expected)
+
+        source = 'tom + dick'
+        expected = ['tom', 'dick']
+        test_it(source, expected)
+
+        source = 'tom * (dick + harry)'
+        expected = ['tom', 'dick', 'harry']
+        test_it(source, expected)
+
+        source = 'tom + dick**0.5 / (harry - tom)'
+        expected = ['tom', 'dick', 'harry']
+        test_it(source, expected)
+
 #-------------------------------------------------------------
 if __name__ == "__main__":